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foreword 


of  related  papers  covering  various 


to  investigate  the  flow- field  variables 


This  report  is  one  of  a  series 

aspects  of  a  broad  program 
associated  with  hypersonic-velocity  projectiles  in  tree  «ight  under 
controlled  environmental  conditions.  The  experimental  research  is 
being  conducted  in  the  Flight  Physics  Range  ot  GM  Defense  Research 

Laboratories,  General  Motor.  Corporation,  and  is  supported  by  the 

under  Contract  No.  DA-01-021- 


Advanced  Research  Projects  Agency 
AMC-H359(Z).  it  ia  intended  that  this  series  of  .eports. 
completed,  will  provide  a  background  of  knowledge  of  the  phenomena 

and  thus  aid  in  a  better  understanding  of 


involved  in  the  basic  study 


the  data  obtained  in  the  investigation. 


ABSTRACT 


A  Fortran  IV  computer  program  for  calculating  nonequUibrium 
streamiube  flows  has  been  developed.  Inviscid  gas  properties  can  be 

s 

obtained  for  quasi-one -dimensional  flows  following  varying  conditions 
of  pressure,  density,  velocity,  or  cross-sectional  area  arbitrarily 
specified  in  the  streamwise  direction.  In  this  program,  the  kinetic 
model  can  include  up  to  20  species  and  40  reveroible  chemical  reactions. 
The  numerical  integration  method  has  a  mathematical  accuracy  check  to 
determine  step  size.  This  eliminates  the  need  for  repeated  trial  runs  for 
optimization.  Controlled  injection  of  small  perturbations  is  employed  to 
run  the  program  with  some  chemical  reactions  near -equilibrium. 
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LIST  OF  MAT  HEM  A  tlCAL  SYMBOLS 


NOTE;  PRIMFS  DENOTE  DIMENSIONAL  QUANTITIES 


4'  4' 

HFi  > 

Ki 


cross-sectional  area  of  stream  tube  «s  £ //l.o 

denotes  vibrational  coupling  of  J  ^  species  to  i  reaction 
boundary  condition  for  the  streamfcube 
number  of  elements 

th  „ 

specific  heat  of  J  species  ”  pT“ 

n  th  •  th 

energy  of  X  electronic  level  of  j  species 

th  ,  th 

energy  of  V~  vibrational  level  of  j  species 

•  th 

change  in  standard  free  energy  for  the  t  reaction 

a  th  -th 

degeneracy  of  ■<  electronic  level  of  j  species 

enthalpy  of  mixture,  //  •=.  •— ~~- 

Ao  To  / 

-  th  .  hj 

enthalpy  of  j  species  including  heat  of  formation  hj  - 

„,f  Rq  T0 

•  th  /  ^ 

energy  of  formation  of  /  species;  h  ■  » — ", 

J  RjTu 

forward  and  backward  reaction  rate  coefficients 


equilibrium  constant  of  (,  reaction 


Boltzmann  constant 


Mi 

MU' ' 


reference  length  in  stream  tube 
*  th 

denotes  j  species 

i  A.  f  f 

molecular  weight  of  gas  mixture;  Mw  s  — — — 

*  th 

number  of  electronic  levels  included  for  j  species 

.  th 

number  of  vibrational  energy  levels  of  j  species 
number  of  atoms  per  species 


Avagadro  number 


f'  p  =  ^(u,f 

C'ij  mole -volumetric  rate  of  production  of  Mj  from  reaction  C  i 

Q;:  =  Qj;  tIgLjL 
»  f/  u;  />; 

$(f)  vibrational  partition  function 

&o  universal  gas  constant  (  Ro  -  1.9B647  cal/mole-°K) 

r  number  of  chemical  reactions  considered 

S  number  of  species  in  the  gas  mixtur 


T'crTj 


translational  temperature;  ~j ~  «=  _ZL_ 

•  th  7"f 

parameter  of  j  species  having  units  of  temperature;  7p •  « — 5- 

»  tli  ^  ^  f 

parameter  of  ^  species  having  units  of  temperature;  "//»*  *  X. 

i,L 

characteristic  probability  parameter^of  species  with 

units  of  temperature;  .  ,  _  Uj 

Uj  »  ~^7 

i  •  .  /  ' 

velocity;  (J  «  — — - 

u: 

,  _  *  4*|i 

vibrational  coupling  factor  of  j  species 

coordinate  along  streamtube:  £/  » 

atoms  of  A.  element  per  molecule  of  j  ^  species 

ratio  of  specific  heats 
4  th 

concentrati  >n  of  j  species  in  moles  per  unit  mass;  *  T^C/WW/) 

th  f  »  « 

vibrational  energy  of  /  spjcies;  6/  «■  — r~, 

J  % 

characteristic  rotational  temperature  of  ^  species; 

"7»/ 

characteristic  vibrational  temperature  of^y  ^  species;  Qvj 

v  T» 

nondimens  ionalizing  term, 


_A_  * 


vibrational  relaxation  distance  of 


V 


*  th  .  4 j 

,  species;  X  *  ■=.-£*- 

J  l  L* 


/{ 


Hi 


chemical  potential  of  j  species 


•#  '  tli  ,  *  th 

stoichiometric  coefficients  of  j  specie s  in  t  reaction 


P' 


^  density;  &  —  — ” 

"t!y  vibrational  relaxation  time  of  j  “*  species 


th 


cy*  .  .  >  th 

degree  of  o.onequilibrium  of  L  reaction 


.  th 


COJ  vibratioh.il  frequency  of  j  molecule 


Subscripts 

GO  refers  to  vibrational  equilibrium 

^  reference  condition 

'  •  .  ,  j  th 

t  pertaining  to  L  reaction 

.  .  th 

j  pertaining  to  j  species 

X  pertaining  to  £  ^  electronic  level 

th 

V~  pertaining  to  is-  vibrational  level 
Primes  denote  dimensional  quantities 


NOTE:  INITIAL  CONDITIONS  MUST  OBEY  EQUATION  OF  STATE  WHICH 
EMPLOYS  THE  CONSTANTS  USED  IN  THE  PROGRAM 


P  T  ^  v 

7"  7-  C  8-  5/  JVo5  x(o  ) 

m. 


I  INTRODUCTION 


Current  studies  of  hypersonic  aerothermodynamics  retire 

rapid  solutions  for  high -temperature  properties  of  gas  mixtures  in 

the  chemical  relaxation  region  of  a  flow  field.  Normal  and  bow  shock 

wave  computer  programs  have  been  written  by  the  Cornell  Aeronautical 
1  2 

Laboratories  '  to  calculate  the  inviscid  nonequilibrium-flov/' field  gas 

properties  behind  a  normal  shock  wave  expanding  alonr  a  constant-area 

stream  tube  or  in  the  shock  layer  of  a  blunt  body  at  hypersonic  speed, 

but  frequently  the  nonequilibrium  gas  properties  must  be  determined 

behind  a  shock  wave  along  an  expanding  stream  tube  with  given  streamwise 

cross-section  area  or  pressure  variation.  Hypersonic  blunt-body  flow 

fields  and  nozzle  flows  of  reacting  gas  are  examples  of  technically 

significant  problems  of  this  nature.  To  fulfill  these  requirements,  an 

irviscid  nonequilibrium  streamtube  flow  computer  program  has  been 

2 

generated  from  CAL'b  normal  shock  wave  computer  program. 

The  frozen  normal  shock  calculations  to  get  the  post- shock 
properties  in  the  CAL  Normal  Shock  Wave  Program  has  been  replaced 
by  reading  in  all  the  initial  conditions  at  the  starting  point  of  the  streamtube 
expansion.  We  allow  an  additional  variable,  the  cross-sectional  area  of  the 
streamtube  denoted  by  the  symbol  A.  This  makes  the  gas  density,  ^  , 

in  the  continuity  equation  a  function  of  both  flow  velocity,  U  ,  and  cross- 
sectional  area,  A.  Hence,  we  need  one  me  e  differently  1  equation,  fee 


equation  of  boundar  y  condition,  to  solve  for  this  additional  variable,  A 
The  boundary  condition  for  this  streamtube  computer  program  can  be 
the  streamwise  variation  of  pressure,  cross-sectional  area,  density 
or  ’velocity. 

The  integration  method,  Runge-Kutta  method,  which  it  used  in  the 
CAL  Normal  Shock  Program,  of  solving  the  ordinary  differential  equations 
i?  modified  with  the  Richardson's  extrapolation  to  increase  the  computing 
speed  (Section  III-A).  With  this  modification  the  running  time  is  reduced 
by  more -t^ian  one  half.  In  addition,  an  accuracy  check  is  included  in  the 
modified  Runge-Kutta  Method.  Therefore,  experience  with  specific 

9 

previous  problems  is  unnecessary  for  the  choice  of  step  size  needed  for 
accurate  results.  This  saves  computing  time  which  would  otherwise 
be  required  for  test  runs  to  judge  the  accuracy  of  the  results. 

The  running  time  for  a  chemical  nonequilibrium-fow -field  calculation 
is  usually  quite  long  in  a  region  where  one  or  more  of  the  reactions  are 
near  chemical  equilibrium.  This  occurs  because  the  concentrations  of 
the  species  which  are  near  chemical  equilibrium  begin  oscillations  around 
the  local  equilibrium  concentration  values  of  that  species;  hence  the 
numerical  integration  step  size  should  be  very  small  to  get  accurate 
results.  An  artificial  perturbation  on  the  degree  of  nonequilibrium,  X*  i  , 
of  each  reaction  is  used  in  this  program  to  save  computing  time  near 
equilibrium.  The  procedure  ia  described  in  Section  III-B. 
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be  used  to  calculate 


This  streamhibe  computer  program  can 
the  nonequilibrium  inviscid  high-tempcrature  properties  of  gas  mixtures 
for  quasi-one -dimensional  Dows  and  car.  be  run  near  equilibrium.  The 
boundary  conditions  of  this  streamtube  computer  program  are  the 
streamwise  variations  of  pressure,  density,  velocity  or  cross-sectional 

area  of  the  flow  field. 


Thi,  report  briefly  discuescs  the  governing  equation,  and  the  numerical 
procedure.,  used  in  this  program.  The  general  operating  procedures  are 
also  presented. 


H.  GOViRfsiNG  EQUATION 


In  the  streamtube 


program,  he  variables  are  velocity  (U).  preeeure 


m’  den8H>'  '  r  >■  cross-sectional  area  (A),  species 
and  vibrational  energy  (  ^  ).  The  variables  U,  p, 
no ndime n s io nal  forms; 


concentration  (  f-  } 
p  ,  and  A  are  all  in 


{2-1} 


(2-2) 


(2-3) 


(2-4) 


PH-ned  qualities  are  dimeneional  and  eubacrip,  denote.  re£erence  conditio, 


The  conservation  equations 


on  a  body-fixed  coordinate 


system  are 


1,2 


Continuity  Equation 

de. _ ?.  c/(J  (2“5> 

^  U  </y  A  tjy 


Momentum  Equ^tiop 


Energy  Equation 


i..  9 


■4ij  iiM ^  Jp 


$  ,  5  S  k*  v  v, 

x  Z %<Cf*~X Z  C  -A  p  t  ^  ^  > 

j(t|  ate)  j=| 


Conservation  of  j  Species 


dlj  «  -L  V  Q. 


j  *  cti/  •- , 


Gons©r>&ti on  of  k  Element 


Boundary  Condition 

The  equation  for  the  streamti'be  boundary  condition  can  be  any  one  of 
the  following: 


dA 

dJL , 


c<8ctcy) 

aBCi(y) 


d  6c3Cy> 


(2-10) 


(2-11) 


(2-12) 


d  u  _  ^  sc4cy) 
73  <*j 


(2-13) 


Where  5Q,(y)  represents  or  A  as  a  function  of  y  either 


in  a  polynomial  or  tabular  form. 


The  temperature  is  computed  from  the  equation  of  state 


p  *  Mw 


(2-14) 


The  reversible  chemical  reactions  in  the  stream  tube  program  involving 
the  specius  Mj  (j  =  1»  2,  3,  . , . ,  5, )  are  represented  by 


$  s  ¥■  w 

vet  %  £  ^ 


(2-15) 


where  an<j  z/Jj  are  the  stoichiometric  coefficients. 


P'j =  ^ 


(2-16) 


P'  "  P'j 


(2-17) 


(2-18) 


The  equilibrium  constant  can  be  expressed  as  a  function  of  the  change 


in  free  energy 


R*Tj 

K :  ~  d?  (&o  To  ) 


(2-19) 


where 


rC  5  £,.° 

~<£~  =  T  n  •;  &- 
/V  rtf  r 


(2-20) 


yi*  j  is  the  chemical  potential  of  the  j  species  at  standard  pressure. 


where 


Based  on  tha  molecular  model  of  a  simple  harmonic  oscillator,  (4j  can 


be  given  as 


4 


fiJit 


where 


>jI _ lUA 

q.  JJ  T 

djl- ! 


i  -  ,-«>■« 


1  4  P-21) 


oj  -  bj  *  r/ 


(2-22) 


6j  M-ty-Ot-e#  *up., 


(2-23) 


In  the  equation  of  conservation  of  j  species,  the  molar  volumetric 

2 

rate  of  production  of  species  j  from  reaction  l  ,  ,  can  be  calculated 

%  ‘  »<  f%  ^  ^  ^ + * 

+*zfrh*  4  rjl]  i 

where  Wx  t  i>x  ,  and  Z2  may  be  1  or  0  and  indicates  the  degree  of 

nor.equiiibrium.  The  first  term  in  Eq.  {2-24)  ( l\]z  =  f /J)J  =>  Bx -0  )  provides 
for  a  particular  nonreactive  collider  in  the  chemical  reactions,  such  as 

Ci  +  O  ^  ZO-hO 

The  second  term  alone  in  Eq.  (2-24)  (Pj  »/,  K/rc2xeC) combines  all  collision 
partners,  such  as  v*ZO+M  .  The  third  term  alone  in  Eq.  (2-24) 

(l  \fij-  =J)  ■=■  O  )  is  included  if  some  subgroup  of  species  are  lumped 
together  as  a  ningle  collision  partner,  such  as  hJO-f"  M  ^  A/  +  O  +  M . 


The  variation  of  vibrational  energy  can  be  expressed  as  follows 


Nonpreferentiai  Model 
Vibrational  Energy  Equation 


+  (/ _ _ f-£  (i 

A  py/rr-y^  t  ~»AteT  f  j  VJ 


a 

r 


7  -  ^  -/? 


(2-25) 


Vibrational  Coupling  r  actor 


V  =±  he - :  <? 

v  Mj  e^/r-j-l 


rt 


-/ 


f 


(2-26) 


B.  Preferential  Model 


Vibrational  Energy  Equation 


W  Xj  \EJ  $]£,  X; 

T  i-Kt 

r;pu  ~x~ 


-fc-sf 


(2-27) 


where 


/ 


£/  *  -  Z.  hi- 

7  &(rFp  v  V  e 

7 


'2-28) 


and  Q.  is  the  vibrational  partition  function. 


Vibrational  Coupling  Faetor 


y  _  Q(r)a(TPi) 

J  “  QXTv)Q(-Uj) 


(2-30) 


The  enthalpy  and  specific  heat  are  given  as 


-r* yi/r) 


h  .fS^C^Q  I  Jt  j  f  ggf. _ * 

J  l  7^  J  '  (?'Vir_  f  &  si  _ 


t  & 


« 

J  =  1,2 . ^  ^  +  **  "  •  /  ^ 

where  V)j  is  1  or  2  and  denotes  the  number  of  atoms  per  molecules. 


6'j/t 


For  vibration  nequilibrium  the  second  tern,  in  Eq.  (2-31)  is 


written,  an 


0v' 

(y\  -I)  — --X.— 

^  VW'. 


and  the  second  term  in  Eq.  (2-32)should  be  omitted. 


■■■■■I 


m  NUMERICAL  PROCEDURE 


The  procedure  for  solving  a  set  of  simultaneous  ordinary  differential 
equations  of  first  order  is  to  solve  for  the  derivatives  of  every  variable 
with  respect  to  y  first  and  then  integrate  streamwise  to  get  the 
properties  at  every  point.  The  modified  method  of  elimination  (Guess’ 
method)  is  used  to  solve  the  set  of  simultaneous  linear  equations 
for  the  derivatives  of  every  unknown  with  respect  toy  .  Then  a  fourth 
order  Runge-Kutta  method  with  Richardson's  modification  (See  Section  III-A) 
is  used  to  integrate  along  the  streamtube.  This  integration  method, 
Runge-Kutta -Richard  a  on' s  Method,  has  the  folldwing  two  advantages: 


1.  It  is  very  easy  to  put  in  a  mathematical  test  for  an  accuracy 
check  without  spending  excessive  computing  time  .  With 
this  check  to  determine  the  starting  integration  step  size  , 
uncertainties  in  the  choice  of  a  proper  initial  integration  step 
size  are  eliminated.  This  saves  considerable  computing  time. 

2.  With  the  residue  modification,  the  integration  step  size  can 
be  much  larger  to  get  results  of  a  given  accuracy.  This  will 
usually  increase  computing  speed  by  a  factor  of  3  or  4. 


In  the  relaxation  zone  of  a  complex  reacting  flow  field  such  as  an 
air  flow,  some  of  the  chemical  reactions  approach  local  quasi-equilibria 
(i,  e.,P£;  0  )•  In  this  region,  the  concentration  of  some  or  all  species 

may  fluctuate  around  the  quasi -equilibrium  concentration  of  that  species. 
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Hence,  the  integration  step  must  be  restricted  to  a  very  small  sise  to 
guarantee  the  accuracy.  Thie  slows  down  the  computation  very  much  in 
the  loc  1  quasi-equilibrium  region.  A  method  involving  artificial  perturbation# 
on  the  chemical  reaction  near  equilibrium  (See  Section  ill-B)  ha#  been 
developed.  With  this  method  the  nonequilibrium  program  is  used  for 
rapid  calculation  through  the  local  quasi-equilibrium  region.  In  the 
perturbation  method,  the  outer  bandwidth  around  zero  uf  _+  0.  1  and  the  inner 
bandwidth  around  zero  of  +  0.  05  are  recommended.  Section  III-B  gives 
definitions  of  the  bandwidths  in  the  description  of  the  method. 

The  boundary  condition  for  this  streamtube  program  is  either  in 
polynomial  or  in  tabular  form.  When  it  is  in  a  tabular  form,  the  derivative 
of  the  boundary  condition  function  is  first  found  with  respect  to  y  by 
a  divided  differences  method.  Then  it  is  integrated  to  get  the  boundary 
condition  function  for  each  y  .  To  solve  the  set  of  simultaneous  differential 
equations,  Eq.  (2-5)  to  Eq.  (2-9),  the  derivative  boundary  condition  function 
with  respect  to  y  must  be  evaluated.  And  also,  the  derivatives  of  every 
variable  are  solved  with  respect  to  y  first.  Integration  is  then  carried 
out  to  get  the  variable  values  for  each  y  .  That  is  why  the  derivative  of 
the  boundary  condition  function  is  obtained  from  the  table  instead  of  inter¬ 
polating  directly  from  the  table  fc-*  each  interval.  In  this  way,  the  boundary 
values  obtained  by  differentiation  and  then  integration  will  usually  not  be 
exactly  the  same  as  the  boundary  condition  values  directly  interpolated 
from  the  table.  For  some  cases,  a  boundary  condition  function,  such  as 


pressure,  may  decrease  with  J  in  a  few  orders  of  magnitude  from  the 
beginning  to  the  end  of  a  region  we  considered.  In  this  case,  the  cumu¬ 
lative  errors  will  be  very  large  near  the  end  point,  if  there  are  very 
small  errors  in  the  beginning  due  to  the  indirect  way  to  get  the  boundary 
condition  value.  To  correct  this  kind  of  error,  the  boundary  condition 
function  is  also  interpolated  from  the  table  directly  at  every  point  and 
is  used  as  the  current  value  instead  of  the  value  from  integration  after 
each  interval.  With  this  correction,  the  boundary  condition  will  follow 
the  prescribed  boundary  condition  function. 


A.  RUNGS 


_kuTT  a- RICHARD  SON  METHOD 


1 .  introduction 

The  Runge-Kutta  method  is  cue  of  the  popular  method,  for  aolviug 
first-order  ordinary  differential  equations  numerically.  t  or  higher- 
or(ier  differentia,  equations,  every  higher-order  differential  equation  i. 
usually  represented  by  a  set  of  first-order  differential  equations.  Thus, 
the  Runge-Kutta  technique  is  also  applicable  to  high-  .order  o 
differentia’  equations. 


When  the  Runge-Kutta  method  is  employed  «  solve  differential 
equations,  the  accuracy  of  the  results  will  vary  inversely  with  the  integration 
Step  sire.  To  ascertain  the  use  of  the  proper  integration  atep  .fee,  after 
each  s^P  is  completed  the  step  sire  is  halved  and  the  integration  over  the 
whole  interval  just  completed  is  repeated.  Then,  by  comparing  the  re.uit 
of  the  full -integration  step  sir.  withthe  one  with  half-step  .Ue  and  integrated 
twice,  we  can  find  the  deviation  of  Ue  result  by  Runge-Kutta  Method  from 
Ue  exact  value  for  this  step  and  see  if  the  step  sire  we  just  used  is  proper 
or  not  in  this  region.  However,  this  hind  of  mathematical  check  of  Ue 
Runge-Kutta  method  usually  triple,  the  computing  time. 


to  aerospace  fields,  many  problem,  require  Ue  numerical  ro.utieu 
of  a  set  of  differential  equations;  for  -ample,  th  calculation  of  Ue  gas 


properties  for  a  nonequUibrium  flow  field  around  a  blunt  body.  For  Ut. 


problem,  CAI,  ha.  a  computer  program.  2  In  Uat  program.  Uey  use  U. 


Mfe  can  use  the  Taylor  series  to  expand  Eq,  (3-2)  as  a  function  si  ^ 


and  X  to  the  N  '  order. 


,  H  rut  _  » 

y®  FCx.J+h  FCx0)  ■+  ***  +  -p~  F  '(X0) 


(3-3) 


where 


X*  =  the  initial  value  o£  X 
F  =  X  ~ 

R-oSF#“”«) 

^  =  a  point  between  X  and  X0 

When  we  use  the  fourth  order  Runge-Kutta  method,  N  *4,  the  remainder,  R, 


we  neglected  is 


The  deviation  from  the  exact  value  for  each  step  is 


(3-4) 


(w).  -  £*£!? 


(3-5) 


if  toe  step  size,  K  ,  we  used  is  not  very  large  we  can  consider  this 
deviation  tc  be  proportional  to  toe  fifth  power  of  the  step  size,  H  ,  only. 


Or  we  can  write 


(Dev)*  **  K  *  h? 


(3-6) 


where  K  is  a  constant. 


If  we  choose  half  the  step  size  and  integrate  it  twise 

Ct>ev)j.«K(Tfx^  *  K^TT 


(3-7) 


1 


TU^  residue,  the  difference  between  the  deviation  from  full 


-grsUon  step  size,  (Dev)^p  to  the  one  with  half  integration  step 


size,  (Dev)  ,  is 

w 

Residue  =(Dev^  -(Dev^ 

=  K(hJ-|; 

°  If  ^  1,5 

-  15 

Hence  subtracting  one  fifteenth  the  residue  from  the  result 
obtained  by  using  half  integration  step  size  and  integrating  twice  will 


(3-8) 


recover  the  exact  value,  if  the  deviation  for  each  step  is  proportional  to 
the  fifth  power  of  the  step  size.  This  modification  is  the  Richardson's 

4 

extrapolation. 


The  procedure  of  using  Richardson's  extrapolation  is  as  follows: 

Suppose  we  know  the  function,  )»  in  Eq.  (3-1),  Using  the 

Runge-Kutta  method,  we  integrate  Eq.  (3-1)  from  XQ  to  X0+h  with  a  full 
integration  step  size,  f")  .  We  get  a  result,  y^y^y^  at  .  If  we 

integrate  Eq,  (3-1)  again  from  XQ  to  %e4  h  with  a  half  integration  step  size  , 
h  A  -  We  may  get  another  result,  +  at  *>X0+h  .  If  the  step 

size,  h  ,  is  small  enough,  y(  may  be  equal  toy^.  But,  usually  they  are  not  equal. 


From  the  definition  of  residue,  Eq.  (3-8),  we  Haw 


Residue,  ^ 

Hence,  with  the  Richardson's  extrapolation,  the  result  should  be 

y  =  y«*  +  Yz  ~(Dev)A 

=  yo  +  y2  -  Residue/ 15 

=  y»  +  &  /is  (3. 

This  result  should  be  very  close  to  the  exact  solution  cf  Eq.  (3  .1). 
When  we  use  the  Runge-Kutta -Richardson  method  to  integrate  a  differential 
equation,  if  we  can  make  Residue/  y,  =  (  y,  -  ^  )/y  <  £  ,  then  th«! 

deviation  of  the  result  from  the  exact  solution  should  be  much  smaller 
than  6  .  Hence,  we  can  use  the  residue  test.  Residue/  V  <6  as  a 

procedure  to  control  the  integration  step  size.  The  residue  teat  takes 
very  little  computing  time. 


3.  Example 


Suppose  we  have  a  differential  equation 

with  the  initial  condition 


y  =  1,  at  X  =  0 

The  solution  of  Eq.  (3-11)  is 

y  =  i+x  -  +  x3 -  X4,  x‘  ,3. 

We  can  use  the  Runge-Kutta  method  and  also  the  one  with  the 
Richards  on*  a  extrapolation  to  eolve  Eq.  (3-11)  numerically  at  a  different 

step  «ire  to  see  how  much  the  Richardeou'a  extrapolation  wiil  improve 
the  results. 


In  Table  1,  AX  =  0.  8,  the  results  obtained  by  strictly  Runge-Kutta 
method  are  accurate  to  the  3rd  digit*  The  results  with  the  Richardsons' 
extrapolation  are  accurate  up  to  6  or  7  digits. 

In  Table  AX  =3.  2,  the  results  obtained  by  the  Runge-Xutta 
method  with  Richardson's  extrapolation  are  even  better  than  the  results 
by  strictly  Runge-Kutta  method  with  AX  =  0.  S. 

4.  Conclusion 

From  the  above  example,  we  can  see  that  in  Table  2,  the  AX 
we  used  is  four  times  the  AX  we  used  in  Table  i.  That  means  the 
Runge-Kutta  method  with  the  Richardsons'  extrapolation  we  used  in 
^abie  2  consumes  the  same  order  of  computing  time  as  the  results  we 
obtained  from  the  unmodified  Runge-Kutta  method  in  Table  1.  Hence,  with 
the  Richardsons'  extrapolation  we  improve  the  results  and  impose  a 
mathematical  check  to  assure  accuracy  with  the  same  order  of  computing 
time. 


The  ‘jpa.ce  needed  for  this  modification  with  the  residue  test  is  very 


small  in  an  actual  computer  program. 


TABLE  1 »  NUMERICAL  SOLUTION  Os^  EG*  <3—1  i  I 


(DX«0*Sl 


X 

O.OOOOOOOE  00 

а. oooooooE-oi 

l ,  f 900000E  00 
2*4000000E  00 
3.2000000E  00 
4.0000000E  00 
4*8000000E  00 
S.6000000E  00 

б. 40C0000E  00 
7 • 2000000E  00 
3.0000090E  00 

e.eooooooE  oo 

9.6000000E  00 


TABLE  2* 

X 

0*00000G0E  00 
3.2000000E  00 
6.40Q0000E  00 
9,60000Q0E  00 


YCEXACT) 
l.OOOOOOOE  00 
K.3279360E  00 
»e*709055SE  00 
-1.3319032E  02 
“89 1 632700E  02 
-3.2749997E  03 
-1.0120038E  04 
—2*6 l 662502  04 
-5,9431 !83E  04 
-1.2232264E  05 
-2.3301495E  05 
—4 • 1 701 388E  05 
-7,0891 059E  05 


NUMERICAL  solut 

Y(EXACT) 

S ,0000000E  00 
-8, i 632700E  02 
-5,9431 183E  04 
-7,08910592  05 


Y (R-K  > 

l.OOOOOOOE  00 
1 ,308821 3E  00 
-0.&128207E  01 
-1 .3344423E  02 
-6*1 679&88E  02 
-3,2757505£  03 
-1,0121 1 37E  04 
-2,61 67761 E  04 
-5.9433173E  04 
-1 o223251 BE  05 
-2,3301 809E  05 
—4  * 1 70 1 773E  05 
-7,089 1 51 8E  05 


ION  OE  EQ*  (3-1 1 ) 

Y(R-K) 

l.OOOOOOOE  00 
-9.3656376E  02 
-5,99401 01 E  04 
-7.100767TE  05 


l.OOOOOOOE  00 
1.3279360E  01 
-8.7090547E  01 
-1.3319032E  02 
-8, 1 63269SE  02 
-3.2749995E  03 
—1,01 20039E  04 
-2.6166251E  04 
-5,9431 184E  04 
-1.2232264E  05 
-2.3301495E  05 
—4 , 1 701 390E  05 
-7,0891 060E  05 


{ DX=3,2 ) 

Y<R”K-R) 
l.OOOOOOOE  00 
-8.1632703E  02 
-5.9431 188E  04 
-7,0891 071 E  05 
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n  nr'nmr?T>t>  jt?rpun^  AM  AlUTMtA  AT 

JL>*  i  .u^v  i,  u^Vu  ilAl  n  Wwtf  A**V>/  wit 

REACTIONS  NEAR  EQUILIBRIUM 


In  a  gas-mixture  flow  field  with  variation  of  thermodynamic  properties, 
the  chemical  reactions  tend  to  maintain  the  gas  mixture  in  chemical 
equilibrium,  Let  us  represent  the  set  of  chemical  reactions  involving  the 
species  Mj  (j  =  ,1, 2,  --- s )  by 


i>.,  M.  5=?  2 


(3-13) 


where  i  =  1,2,3  ---  r  is  the  number  ot  reactions.  V;  .  andj*-*  are 

AJ  AJ 

f  f 

stoichiometric  coefficient?  for  reaction  i.  ^  and  .  are  forward 

and  reverse  rate  constant  respectively.  Tht  >4p^.  can  be  calculated 

from  the  local  vibrational  equilibrium  rate,  , 

>ia> 


v  (k)a,j 

r+  '•a,  j- 1  J 


(3-i4) 


where  '/j  is  the  vibrational  coupling  factor  and  may  be  1  or  0  ,  denoting 
which  reaction  I.  will  be  affected  by  the  vibration-dissociation  coupling  process 

The  net  volumetric  rate  of  production  of  species  i  from  reaction  i  <?!". 

'  \) ; 

is  given  by 


<Sf4  -■($  [if  f*  7T  1^-  ^ '  o  A  77  r  /  n^L 

Aaj  (v  a  Jo,*-* 


P-15) 
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wners 


Pi  *  <-r.  P<u  0*o f  Pi  *  p. 

gi  -  ■  at  “/ 


# 


Since  the  equilibrium  constant 


Eq.(3-15)may  be  rewritten 


{  3  17) 


and 


{ 3  ~  1  &) 


The  X-A  is  referred  to  as  the  degree  of  nonequilibrium  of  the  ^ 
reaction.  Wlien  the  reaction  i  is  in  equilibrium  *0.  Otherwise,  the  • 
will  have  the  values  between  1  and  -  co,  excluding  0. 


In  the  flow  field  of  gas  mixture,  if  the  chemical  reaction  time  is  slower  than 
the  time  characterizing  the  variation  of  thermodynamic  properties  in  the  flow, 
the  concentrations  of  the  species  near  equilibrium  usually  oscillate  around 
the  equilibrium  values  as  shown  in  Figure  1, 


Figure  1  Variation  of  Species  Concentration  Near  Equilibrium 

Using  the  Runge-Kutta  method  to  calculate  the  concentratiai 

of  species  along  the  streamwise  distance,  y  ,  ihe  integration  step  size 

has  to  be  smaller  than  the  distance  between  two  adjacent  relative 

maximum  or  relative  minimum  points  in  the  actual  concentration  curve 

to  attain  the  desired  accuracy.  Near  equilibrium,  the  distance  in  the  y- 

direction  between  two  relative  extremes  of  the  species  concentration 

curve  is  usually  very  small  and  derivatives  become  large.  Hence, 

the  computer  running  time  for  a  flow  field  at  near -equilibrium  is  very  long. 
ry  , 

The  t  -  variation  of  those  reactions  near  equilibrium  is  very  slow 
along  the  streamwise  distance,  y  .  Also,  in  this  region  the  values  of 
are  very  near  zero.  Ae  an  approximation,  we  can  set  equal 
to  zero.  Since  the  XK  vs  j  curve  is  nearly  flat,  the  results  of  this 
approximation  will  not  deviate  too  much  from  the  actual  values,  but  with 
this  approximation  the  integration  step  size  can  be  ■nuch  larger.  The 
scheme  of  ihis  perturbation  method  is  as  follows. 


Two  specified  bandwidth*.  2 $  and  2€  are  placed  around  sCi  s0 

/V  c  ^  . 

the  Aj't  vs  y  space  whet?  <6  is  less  than  fe  ,  We  ’'.eep  checking  the 

%i  -values  at  every  integration  interval.  If  the  previous  %-i  -  value  is 

not  aero  and  the  current.  !&**.  I  calculated  from  itff  algebraic  definition 

(see  Eq.  (3-18))  is  smaller  than  S  ,  we  set  =0  which  makes  the  $  (J  =0 

from  Eq.  (3-17).  This  term  is  omitted  from  the  summation  in  Eq.  (2-8)  for 

integration,  and  the  reaction  it  represents  is  temporarily  and  artificially 

frozen.  If  the  IXJ  calculated  from  Eq.  (3-18)  is  equal  to  or  greater  than  &  , 

then  the  currently  integrated  %  1  -value  is  used.  If  the  previous  -value 

is  zero,  we  Bet  %[=0  for  IX;|  less  than  £  .  Otherwise,  the  calculated 

value  from  Eq.  (3-18)  is  used  as  the  updating  %l  -value.  Thus,  those  nonzero 

1^-tl  -values  less  than  €  never  enter  the  integration.  When  the  chemical 

reaction  moves  away  from  equilibrium,  calculated  values  from  Eq.  (3-18) 

are  used  until  its  exceeds  1^1  .  In  this  way,  instead  of  using  one 

$  -band,  we  can  avoid  the  rapid  fluctuation  of  the  -values  from  0  to  the 

calculated  value  greater  than  b  and  vice  versa.  This  X/'L~ test  is  an 

empirical  method.  Kence,  the  bandwidths  for  S  and  €  should  be 

determined  by  experiments  for  each  particular  class  of  problem.  Experience 

is  gained  rapidly  and  fixed  £  or  €  values  work  for  a  wide  range  of 

conditions.  Far  less  high-speed  memory  storage  is  required  for  this  method 

than  for  a  separate  analytical  perturbation  code.  With  the  stre&mtube 

program  described  above,  most  of  the  available  cells  are  occupied  so  that 

the  artificial  perturbation  C  i  -  te  s t  represents  a  significant  advantage.  It 

probably  obviates  the  need  for  a  chain  program  modification  in  this  case. 
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IV  PROGRAM  DESCRIPTION 


This  program  is  written  in  FORTRAN  IV  for  an  IBM  7040-44 
digital  computer.  This  program  requires  Approximate.y  25,ClO  memory 
cells  and  at  least  two  magnetic  tapes  for  the  normal  outputs.  The  system 
at  GM  DRL  is  composed  of  a  4K  central  system  monitor,  with  32K  of 
core  storage,  8  magnetic-tape  drives  and  a  1611  card  read-punch  unit  on 
channel  A.  With  specific  sense  switch  settings,  we  can  select  optional 
outputs  or  tape  #2  and  tape  #3,  The  sense  switch  settings  and  tape 
requirements  are  as  follows: 

\ 

A.  Sense  Switch  Settings:  . 

#1  -  optional  output  on  tape  #2  of  species  thermodynamic 

and  chemical  variables 

#3  -  optional  output  on  tape  #6  of  the  array  of  coefficients  r 
the  set  of  simultaneous  differential  equations 

H  -  optional  output  on  tape  #6  of  the  derivatives  of  boundary 
condition  function  with  respect  to  y  . 

#5  -  optional  output  on  tape  &  of  the  molar  rate  of  production, 

#6  include  %i  (  i,  =  >  •  '  )  in  each  result  paragraph 

on  tape  #4. 

B.  Tape  Requirements: 

Output  tape  #2  (optional  output).  Tape  #2  is  written  only 

when  sense  switch  #1  is  on  and  DUMP^  (input  variable) 
is  set  equal  to  1  or  4,  A  paragraph  of  species  thermodynamic 
and  cticinical  variable  is  printed  for  the  first  Runge-Kutta 
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a$f*T%  Hiryt)  r  \  Of  all  a***-**.* 

~  — nsjx)  J  —  ~~r~, 

DUMP  =4  in  an  interval. 

Output  Tape  #3  {optional  output}  -  Tape  #3  ie  written  only 
when  sense  switch  #5  is  on.  For  each  success  step  a 
paragraph  of  the  molar  rat*  of  production  of  species 
(  j— C S  )  from  all  reactions  t  ,  &ij  ,  is  printed  out. 

Output  tape  #4  (normal  output)  -  A  paragraph  of  results  for 
each  printing  invervai  in  printed  out. 

Output  tape  #6  (normal  output)  -  mirror  messages  of  failing 

test  are  printed  out  on  tape  #6.  Whea  sense  switch  4  is  ou, 
the  derivatives  of  boundary  condition  with  respect  to  y 
are  printed  on  tape  #6. 

The  chemical  system  of  this  atreamtube  program  has  a  capability 

to  read  in  forty  reversible  chemical  reactions  ^  —  |  ^  3  y  j 

and  twenty  chemical  species  (  j  **!,!,  3,  S  ).  The  species  are 

subdivided  into  independent  species,  the  eiements  {  -ft  =■  (^  2. ,  ,  C  )* 

and  dependent  species  (  j  «C+I,  •  -  •  ✓  5  ).  In  the  dependent  species 

# 

we -can  have  some  diatomic  species  {  j  =.  -p-f  |  ,  .  .  y  ^  )  for  which 

vibration  noneqcdlibrium  are  considered. 

With  this  program,  we  c  an  put.  more  than  one  set  of  ftreamtube 
calculations  in  one  run.  There  is  one  data  card  right  after  the  fis\  set 


of  data  we  have  for  this  Jfarj  w@ 


ot  data  cards  to  tell  hs%v  many  ^eis 
twt  ail  ih«  sets  of  rata  right  after  n  in  ore  deck  and  run  them.  When  et 
1 4  finish***-  th-_-  will  put  ,!end  of  file"  on  every  tape  used.  For  every 

stresrriube  there  »■*&  two  normal  termination  conditions.  One  is 
Y  $  (input  variable).  The  other  one  is  ^  Tm.n, 

(h'cput  variable  TEST,).  At  normal  termination,  the  following  line  will  be 

4. 

printed  c-t-  tape  “6  : 


FUJI .  \run  no. ) 


COMPLETED 


There  are  two  possible  error  terminations.  One  is  due  to  the 
zero  d-terminant  encountereo  in  the  set  of  simultaneous  linear  differential 
equations.  The  other  one  is  due  to  Ayn  ^  A'j  m)n  .  If  either 
occurs,  the  error  messages  are  printed  out  on  tape  #6  and  the  system 
will  go  on  to  process  the  next  set  of  data.  If  the  error  terminate  t  occurs 
in  the  final  set  of  data,  this  run  is  finished. 


Paul  V.  Marrone,  "Inviscid,  Nonequilibrium  Flow  Behind  Bow 
and  Normal  Shock  Waves,  Part  I;  General  Analysis  and  Numerical 
Solutions,"  CAL  Rept,  No.  QM-l626-A-12(I)f  May  1963 
Leonard  J.  Garr,  and  Paul  V.  Marrone,  "Inviscid,  Notiequilibxium 
Flow  Behind  Bow  and  Normal  Shock  Waves,  Part  II:  The  IBM  704 
Computer  Programs,"  CAL  Rapt,  No,  QM-1626-A- 12(11),  May  If  3 
P.  Henrici,  Discrete  Variable  Methods  in  Ordinary  Differential 
Equations,  John  Wiley  &  Sons,  Inc. ,  1962 

±*‘.  B.  Hildebrand,  Introduction  to  Numerical  AnalyBiB,McGraw-Hill 


Co.,  1956 


APPENDIX  A 
INPUT  FORMATS 


A. 

First  Five  Cards 

Page 

A-l 

B. 

Species  Thermodynamics 

A- 3 

C. 

Chemical  Reaction  System 

A- 3 

T>t 

0!  jl.  Matrix 

A- 4 

L. 

Initial  Condition 

A- 4 

F. 

Boundary  Condition 

A- 5 

G. 

Boundary  Condition  Pol’  normal 

A- 5 

H. 

Boundary  Condition  Table 

A-6 

I 


Number  of  Data  Sets 


A- 6 


appendix  a  input  formats 


lorBI«  of  N»t  data  cards  are  giv5a  fa  detail 


below. 


A-  First  Five  Card* 

1  Card  “  Format  (AC),  8l3,  2E10.  4) 

IRUN  run  number 

number  of  species  in  the  gas  mixture 

r  n"mbsr  o£  chemical  reaction,  in  the  ga,  mixture 

f' 3  'Uiine  v^brational  nonoquilibrium  group  of  .pecie. 


number  of  elements 


B-  C.  boundary  condition,  UU.iyo  3  p*  4  „  A 

E'  E1NI,  electronic  excitation  of  the  .pecie.,  0-.no,  1  ^yee 

DUMPM)  °!’tioni‘1  “"‘Pet  °<  itt  .tep  of  Runge-Kutta,  0  -.none, 
1  -* first,  4  *  all 

A y  Sw+.  8 Parting  interval,  in  cm 

t 

^s+op  value  of  y  to  stop  computation,  in  cm 

Second  Card  Format  (3E9.  5) 

/ 

c  reference  pressure  in  dyne'cnA 

/ 

(°0  reference  density  in  gm/  m3 

T' 

o  reference  temperature  in  °K 

Blank  available  for  further  information 

A0  initial  cross-sectional  area  of  streamtube  in  cm2 


A 


! 


reference  velocity  in  cm,*  b@s 
reference  molecular  weight,,  in  gm/moie 
reference  length  In  cm,  usually  one  cm 


Third  Card  Format  (5E14.7) 


TEST. 


TEST, 


TEST, 


TEST , 


TEST, 


number  of  output  result  paragraph*  per  page 
value  of  temperature  that  terminates  the  computation 
percentage  change  in  T  allowed  for  each  interval 
halfwidth  of  outer  X'i  test  band,  usually  =  0.  1 
half  width  of  inner  %;test  band,  usually  =  0.  05 


Fourth  Card  Format  (5E14.7) 


TEST, 


TEST, 


blank 

blank 


TESTg 

blank 

TESTg 

first  val  ue  of  y  to  be  printed  after 
nondimensional.y^- /  L 

TEST  IQ 

t  .  f 

printing  interval ,  &  Jp~A y^,  /L 

Fifth  Card 

Format  (5E14.7) 

TESTU 

number  of  successful  steps  before 

TEST  !  ^ 

C,  step  size  increment  factor 

testi3. 

blank 

testh 

minimum  Ay  allowed 

test15 

maximum  Ay  allov.'ed 

B.  Species  Thermodynamics 


There  are  lour  cards  for  each  of  the  j  =  1,  3  epeeies 

First  Card  Format  (6£  12,  6) 

number  of  atoms  per  molecule,  one  or  two 


J 

b. 

J 


f?v. 

J 


N 

J 


Blank 


Second  Card 


constant  for  chemical  potential  calculation 

characteristic  vibrational  temperature,  in  °K 

heat  of  formation,  calories  per  mole 

number  of  vibrational  levels  for  harmonic  oscillator  cut 

off  at  dissociation  energy 

available  for  further  species  description 

Format  (4F12.6,  4A6) 

describe  vibrational  relaxation  time  for  species  in 
the  £  +  1  -*»  g  group; 


■TV  =  T„:  (t‘)  exf‘[ 


To 


(T-fl' 


J 


oi/x*5 

— 1 *■ - —  Jr<r 

Csn 1 


C. 


SPECJK 
Third  C  ird 

'“jit 

Fourth  Card 


j**1  species  identification 
Format  (12,  8F2.0) 

number  of  elccironic  levels  of  j**1  species, 
degeneracy  of  each  l  ^  electronic  level 
Format  {8E9-  5) 


J?=  1,  — ,  8 


energy  levels  of  electronic  states  in  cal/ molecule 
Chemical  Reaction  System 


There  are  two  cards  for  each  i  =  1,  r  reaction. 


First  Card 


Wi,2 1,  J>3 

J/-* 


Format  {.sF’2,  G%  20FI.Q*  26”  l.S,  EOF  iaO* 

denote,  which  &'.  rri*“™  *®  «« 

-J 

right-hand  side  stoichiometric  coefficient  of 
species  j  on  reaction  ^ 

left-hand  side  stoichiometric  coefficient  of  species  j 


Second  Card 
KFI1ND 


on  reaction  £, 

denotes  which  species  is  involved  with  vibrational 
dissociation  coupling  of  reaction^  1  **-ye8,  0-^no 
format  {14,  4E14.7,  2A6) 

denotes  direction  of  input  reaction  rate  constant, 

Q  -*-forward,  1  -^-backward 


Cfi'  D# >o  reaction  rate  constant 

* 


4’  M  e*p[-  -~sr  }~ 

D«t  r-cc  r  ^  {"7 •'/**»]  fndr-Sr 


Cm* 


Sc-  c  °r  tyxJr-  -  5#tf 


.til 


3PECIK  identification  of  i  reaction 
I>.  jk  Matrix  Format  {20F2.  0) 

Need  one  card  for  each  &=  1, - ,  c  element  to  describe  all 

species  by  the  elements 
E.  Initial  Condition 

First  Card  Format  (5E14. 1) 

Ore  card  for  every  five  of  £  +  l  -s*g  species.  If  f  -  g,  omit  this  card 


6} 


.  "h 


vibiational  energy  ofj  1  species 


A- 5 


Second  Card  Format  (5 El 4. 7) 

Dne  card  for  e^ery  five  coefficients 

COP  coefficient  of  polynomial,  starting  with  the  lowest 

degree  coefficient 

Boundary  Condition  Table  (IPOT*!}  Format  (2&14.  7/ 

One  card  for  each  TY 
TY  coordinate  value 

corresponding  value  of  the  boundary  condition  function 
Number  of  Data  Sets  Format  (12) 


1ST  OP 


number  of  data  sets  for  this  run.  Only  need  to 
put  one  "number  of  data  Beta"  card  behind  the  first 
set  of  data. 


APPENDIX  B 


FLOW  DIAGRAM 


SOI 

MAIN  PROGRAM 

Em* 

B-3 

SO  2 

Subroutine  START 

B-6 

SO  3 

Subroutine  REABIN 

R-7 

SO  4 

Subroutine  WRITE  1 

B-fi 

SOS 

Subroutine  WRITE  2 

B-9 

SO  6 

Subroutine  MATRIX 

B-1G 

SO  7  A 

Subroutine  DER 

^-l\ 

S07B 

PREFERENTIAL  FLOW  I,  II,  III 

B- 15 

S08 

Subroutine  SUB1 

B-18 

S09 

Subroutine  SUB2 

E-19 

SIO 

Subroutine  SUB4 

B-20 

Sll 

Subroutine  SUBS 

B*21 

S12 

Subroutine  SUB6 

B-22 

SI  3 

Subroutine  TEST1 

B-23 

S14 

Subroutine  TEST2 

B-24 

S15 

Subroutine  TEST  3 

B-25 

S16 

Subroutine  QVRFLW 

B-26 

S17 

Subroutine  LCOUNT 

B-27 

S18 

Subroutine  STOP 

B-Z8 

SI9 

Subroutine  SXMSOL 

B»29 

S20 

Function  DERT 

B~30 

S21 

Function  YINT 

B-31 

Processing  or  operational  box.  Descriptive 
or  working  block  of  instructions . 


Call  box.  To  call  a  library  or  closed  subroutine 


MAIN  PROGRAM 


$ccc)  (I rid! 


/Q\ 


t5TA!<T  K  P 
ASSIGN  ItoC-5ooo 


5T.  vRT 


I  i.T  +  )  ;  **  4“  4 


READIN 


WRITE  I 


CALL 

MATRIX 


xeXT 

ii  *  testi( 

(C 4  i) 

«=  c.  ■+  1 

C-f+1) 

=  f  +  i 

C3+0 

-  g  +  i 

f 

-  £ 

M  »4-* 

Ml 

M+  1 

MX- 

S-tg-f 

MX  1  - 

MX+  1 

MX  2- 

MX +2 

MX  3* 

MX  +  3 

MX4- 

MX*  4 

MX  5  « 

k  5X  •*  5 

(f.CZZ  t 


A.I^SoWETiRn 
v  “ 


F0R  je  l,J,( 

Qvj  *  e'VJ  /t,' 

-  h]7(RiT0') 

7'  =  b;  +  {5  •  z  (>?•  - 1  ».M%)/z 


{Page  i  of  31 


55go  _ _ 

I  F«B  l  «  1 ,  r,  t 

fS-Lj  -Vij  -Vu  j*ljS,,i 

%  & 

~J;  -  r  }».• 

■  rr, "  -j 


J'!  ' 

j"1 


F0R  j  *r  1  .  J,  | 

Ej«“Eji/(3.Z982E*2^*T;)  **15,1 
EEC  (hj) «  o 
E£C(yjj/T)  =  o 
EEC(Cfj)  *  0 


<  C.-  R.T0'f4. 13228  £  j  ) 


in  tlx  WINT 

P-j  >1,3,1 


X  *  S'/|.0 


&Y  ~  T83T 


WRITE  flSUT  <&N  TAPE  6 
RUN  SKIPPED 


3*9*1 


mm 


■SWWffsriA 


(Page 


tip 


F0R  j=!,$ 


MW  =  MW  U  -  U 
£-£  ?  =  P 

A  =  A 

’Si-Yj  j-W 


J  NC 


l  _  fcsciah  A  S 

I  1  "  ~T~  ! 


vc  » IC  +  I  i-^<^  IC=NC 


F®R  j  -f+1,9,1 


£j  +  P’t-M  n  Vt,  L 


Ftfft  j  =  1,5, 


9-(F+0 


IFAIL^O 


y»y  +  ?aiAy 

Mw  =  i/^Yj 


4;il -ay  2^i  where 
T1  ;  d*  *,=  U 

ZolUiRiu,  #*«P 


:ii=4/ 

>-— -w-  ili  =  iii+  i  j- 

^£5 

feSSO 

^  P0R  .4+1,3, 

>0  4:  =  €;+£€: 

J  J  j 

If.  _  .  _ 

1  -  ?«mw*  A/P 


ViTj 

Uj-TiKVi-T(j)/|5 

U=U  +  (U~U,)/l5 

?=e-He-?,yi5 

Vl  5 


rajL-£=igS-5S=g 


SUBROUTINE  READ  IN 


READ  F«R  3  **  SPECIES  THE  F0v_LsSWIM<t  4  CARDS 

Tlj  ,  bj  ,  tWj  ,  Kj^  Nj  ,  X'«j 

TAjy  fBj;  ,  (24  COLUMNS  <*F  IDENTIFICATKSN) 

9ji  1-I.0J 


j*j*H 


FOR  EACH  RSAC7I0N  l  l  *  I ,  r,  \  READ  2  CARDS 
Wi.Z^Di,^  ]=!^°'l)/(^j  j-  l,20/|)(Aij  j»  1^0,1) 
KHINDi,  A+iee>  Z  W0«DS*F  IDENTIFICATION 


1520 


F0R  t:ACH  -k  -A-I.C,,!  READ  s  CARD 


READ 

INITIAL  C<^DITIC-M5 

A/ 

ej 

j  *  C-f+t  )  ,  9  ,  S  (IP  3  >  'f  T 

'X* 

■»j 

j  *  i  ,  s  ,  i 

■V^ 

F  F  U  ,  MW 

>  '  / 

READ  ft,  C.  P0LY. 
F4>R  I  “  I  /  N^SR,  | 
IF^P^ReSr^RENW 
c#Pr,  j » 


READ  NOR,!!*! 


READ  B.C.ms 
FOR  i  *  i , , 8 


WRITE  ^TpUT  TAPE  *4  ALL  IfciPUT  »>ATA 
“  l>*j  ,  Ajj  IN  FIXED  P0EHT  FORMAT)  _  _ 


READ  FIRST  FIVE  CARDS 

RuW*j  5,  r,  f  ,  0,  c  ,  XB0P  ,  EE  ,  DUMPV#J , ^y,TArr , 
f  0' ,  pc'  >  T0' ,  T0  y  U„'  ,  MWo  ,  l' 

TSST?  ;  =  :  15  I  (5  PER  CARO) 


100 


1 

ENTRY} - — 

—  .  . . . — " 

WRITE  HEADINGS 

0N  0c;tput  tape  *6  j 

!02 

WRITE  HEADIN&S 

AMD  OUTPUT  FORMAT 

(ZSN  OUTPUT  TAPE  *4 

J 


_ 1 

[return 

SUBROUTINE  DEE 


(Page  I  at  4) 


£NTiv?  L 


_PREPERENTIAL 

preferential! 

FLOW  I 


J  CALL 

\  SUB  1 

*&?>Z 

\tc  j* 

f  CALL 
VSU6  w 


K  T 

t7  *  t;  t 


F0R  j  -  l/  ■F/  I 

^  -  -*j- /.S  +  >-*  -  £EC  (£*-) 

7ij«  !  Cft=  ?.5+EEC(Cpj> 

€j  »  0.0 

Vbj  »  2.5T  +  EEC(Hj) 

r,jt, 

.  EXP  (-=■*■) 

cfr  — ^rr~  +  EECfCpj) 

i>j-  Y(Sti;nj.-i))T-.(r,i-!)£j  +h?  +  EEC  (hj) 


1210  1350 


SAME  AS  1050  -u  90 
Except  rz5r  range  $p  j 

J  “9fLS,t  IN5T6AO  $F  j  = ! ,4, 1 


i^SV3 


^  ---V  *  ■£  i-J- ,T 


9Sc--5ca 


■■HI 


SUBROUTINE 


CALL 


(5+2(nj"  l)HE£C  (Cfj) 

\i  -  j (5 + 2Cnj -I))  J»T + (*»j- 1) M  J +  ~  -  &C  (  y) 

<(i-T] 

(5+2(nj  -l))T 4 (Vij  -|)£j+  bj  4  EEC  (bj) 


Tvj  T 


'  UcU_ 
i  L' 


V  rs  I  !--W*540)  i 

ilzIj  i 


T0 


i-sxp(-Niw;)][exp^)-  l 

Nj[6XP(W;}-  l]K%)-l 


i=j+ 1 


TQ  BE  REPLACED 
BY  PREFERENTIAL 
FLOW  E ,  F^R 

PREFERENTIAL 

MODEL 


--l 


!^9sm9S9i9PHIi^p^aHp9PHB9^pipB 


PREFERENTIAL  FLOW  I 


6MTRY 


YPREV 


REAP  A  CARD  WITH  W>LL0WrMS.  VARIABLES 

MJ,  UJ,  WCj,  W*  fcj,  We  faj 

WRITE  7ME  AB#VE  VARIABLES  «N  0UTPUT  TAP*  #4- 


Nj-Nj+I 


KK  *  I 


Ei/KK)  =  (KK-0,5)(Wej  +(KK-  O  SX-uJeXej  +'KK-C.5)( w,^-  +W#^  (KX-0.5)))) 

Evj(KK}«<  EvjCKtO  -E  vj(!) 

TEMP2*  g  EOj(KK)  EYPLEvj(KK)/Uj] 

TEMP3  *  g  BXPCEvj(KK)/Uj] 


KK-KK4I 


.987522  «  TEMP2 


To  R»'fTEHP3) 
« TEMPS 


-■vi  •• 


PKEFERENTLAL  FLOW  II 


TtMP2  *  rW-(KK)EXP  [-Evj(KK)/^;] 

Hi 

TEMP3=L  EXP[-Ev:(KK)/t;.] 

KK»I 

M; 

TEHP4-£  EXP[-Ev;CKK)/t'3 

KK»I  J 

94* 

TEMP5«i:  &XP[-Evj(KK)/(TvJTt,/) 


SUBROUTINE  SUB1 


StTBROUTENE  SUB4 


WRITE  tiorpOT  TAPE  *2 


WRITE  OUTPUT  TAPE  *2 

<2ij 


WRITE  0’JTPUT  TAPE  *2 

THE  DERIVATIVES  0P  THE  F*tL0WIW€r  VARIABLES  • 

€j  j-^+iygj  Cip  it  exists);'?,  A*,4>;p  ;u 


suBRouTmir  su&6 


tt2  -f>e*p[-^] 

1-1 

T  T I  -  B'.p  [-  Eji 

TT3-^9j(EXP[--itt]E/, 


TT3 


EEC(hj)-  -Ili- 

eec©»j>- 


SUBROUTINE  LCOUNT 


ENTRY 


LINE  CSIONT  -  55 


LINE  CjSUMT  **  0 


WRITE  OUTPUT  TAPE  *  4 
RECORD  T0  EJECT  RA&E 


RETURN 


SUBROUTINE  STOP 


appendix  c 

PROGRAM  LISTING  AND  SYMBOLS 


1.  LIST  OF  FORTRAN  SYMBOLS  C-l 

2.  PROGRAM  LISTINGS  C-9 

501  MAIN  PROGRAM  ^“12 

502  Subroutine  START  C-l 6 

50 3  Subroutine  READIN  C-l 7 

504  Subroutine  WRITE  1  G-21 

505  Subroutine  WRITE  2  G-24 

50 6  Subroutine  MATRIX  C-25 

S07A  Subroutine  DER  C-26 

S07B  PREFERENTIAL  FLOW  I,  IIf  IH  C-30 

508  Subroutine  SUBl  C-33 

509  Subroutine  SUB2  C-36 

510  Subroutine  SUB4  C-37 

511  Subroutine  SUBS  C-38 

512  Subroutine  SUB6  C-39 

513  Subroutine  TEST1  C-40 

514  Subroutine  TEST2  C-41 

51 5  Subroutine  TEST3  C-42 

516  Subroutine  OVRFLW  C-43 

517  Subroutine  LCOUNT  C-44 

518  Subroutine  STOP  C-45 

519  Subroutine  SIMSOL  C-46 

520  Function  DERT  C-48 

521  Function  YINT  C-4S 

522  Common  and  Dimension  Statements  C-50 


APPENDIX  C  PROGRAM  LISTING  AND  SYLSSOLS 


List  of  Fortran  Symbols 
The  Fortran  Symbols  listed  hare  are  j;iyen  in  alphabetical  orAer, 
The  corresponding  report  symbols  and  their  meaning  are  also 


given  in  this  section. 


AH  C 3 > 


Alio)*  y 

A  I  IT$,>  *  MW 


initial  condition  for  a  r^k-r  integration  interval 


*1 IC<3> 

A!T 4  3  > 

A I T  |  ( 3  j 

kcl  (4  ) 

A21CU) 

A2T (4 ) 

L2T I <4) 

A3I <U> 

A3ICCJU 

A3TIU) 

C»2 


Alice  i  )«y 
Au.o:/>-  T 
A|ICC3)“M»v 


UPOAT^o  INITIAL  CONDITION  IN  A  R-JC-R  INTERVAL 


At  TU> - y 
A!  TCI)  *  T 
A I  TO)  ~  MiV 


UPDATED  VALUE  IN  A  R-X-R  INTEGRATION  INTERVAL 


Ai  Ti(l)  -  JT| 
AlTlf/5  T, 
AlTI  (3)  *  MW, 


UPDATED  VALUE  IN  A  R~K~R  INTERVAL*  PULL  I  NT#  STEP  SIZE 


A  2  I C I )  »  C 
A2I(2) 

A.2 1  (3)  =  ? 
A2U4)-X 


INITIAL  CONDI  TION  FOR  A  R-K-R  INTEGRATION  INTERVAL 


A/ICC!)  -U 

A2ICC2)  **  p 
A‘2IC(3)«^ 
AZlCf4)-^ 


A 2  Tv  I)  «* U 
A2T(?)-f 
A2TO)** 
A2T(4)  -  A 


A2T|(I)«U, 
A2TI  (2)“p, 
A2TI  (3J«f, 
A2TK4J-A, 


UPDATED  IWiTtAL  CONDITION  IN  A  R-K-ft  INTERVAL 


UPDATED  VALUE  IN  A  R-K-R  INTEGRATION  INTERVAL 


UNDATED  VALUE  IN  A  R-X-R  INTERVAL*  FULL  INT.  STEP  SIZE 


tij  J  * 


INITIAL  CONDI T! CM  FOR  A  R-ff-R  INTEGRATION  INTERVAL 
UPDATED  INITIAL  COHO! TION  IN  A  R-K-R  INTERVAL 
UPDATED  VALUE  IN  A  R-K-R  INTEGRATION  I NTT  VAL 


*  J~  7%A,_ 

' 

ASTI iJ 1 

i 

« 

t  |P  1  A 

i  *an* *m#A  uAi  « *-*  *  *  ±*im»**i\*.m*  «» -  *  **— ■  -i--. — 

v»4.vw  A  H*  **  r  *s$UlU  ff*f«  BlSf* 

AAttJi 

INITIAL  CONDITION  FOR  A  R-K-R  INTEGRATION  INTERVAL 

A4IC <U» 

j 

UPDATED  INITIAL  CONDITION  IN  A  R-K-ft  INTERVAL 

A AT i J ) 

€i 

UPDATED  VALl'E  IN  A  R-K-R  INTEGRATION  INTERVAL 

A4T1 (J» 

4ij 

j  “f+ 

UPDATED  VALUE  IN  A  R-fc-R  INTERVAL*  FULL  I  NT*  ST)3»  SIZE 

A&F  u  n 

H- 

^los 

FORWARD  03  BACKWARD  RATE  CONSTANT  COEFFICIENT 

ALPHA  14  } 

HUNGE-KUTYA >RICHARDSON  COEFFICIENT 

ALPf Jl J«Kf 

NUMBER  OF  Je  ATOMS  IN  j  SPECIES 

S<4S*45> 

ARRAY  OF  CONSTANT  COEFFICIENTS  OF  THE  SET  OF  SlfcUL*  EOS 

SETA  <41 

runse-kutta-richardson  coefficient 

tSETAIf  |  J 

Pi 

S 

Pi  »£  Kj 

j*l  J 

|  *  SETA?  JU  «J) 

ftj "  2'5  ~^j 

BKFf  If) 

FORWARD  OR  BACKWARD  RATE  CONSTANT  COEFFICIENT 

’  C 1 4S  >43  > 

ARRAY  OF  COEFFICIENTS  OF  THE  5ET  OF  SIMULTANEOUS  EOS. 

cao;© 

A' 

initial  cross-sectional  area 

]  C  A I J  C  I  ♦  U ) 

*3 

VIBRATION-DISSOCIATION  COUPLING  INDICATOR 

'  CC  PJUt 

CPJ 

SPECIFIC  MEAT 

cot  f  n 

3i 

DEMOTES  WHICH  SPECIES  PRODUCTION  RELATION  TO  USE 

CEJLPf J*L> 

Ejl 

DIMENSIONLESS  ENERGY  OF  ELECTRONIC  LEVEL  l 

€E„ffLPKf*J*L> 

ENERGY  OF  ELECTRONIC  LEVEL  £ 

■^hi  id) 

oox 

DEGREE  OF  MONEOUILISRIUM  OF  REACTION  t 

CHI  ITS n 

1 

DEGREE  OF  NONECUl LIBRIUM  OF  REACTION  i  ON  PREVIOUS  STEP 

ocp  i  a  i 

i 

FORWARD  OR  BACKWARD  RATE  CONSTANT  COEFFICIENT 

ckkii 

K* 

EQUILIBRIUM  CONSTANT  OF  REACTION  i 

CKP  MU 

«K 

\ 

EQUILIBRIUM  CONSTANT  OF  REACTION  L 

1 

CLP 

L' 

REFERENCE  LENGTH 

eiesep 

WWc 

REFERENCE  MOLECULAR  WEIGHT 

* 

C-3 

l 

fi  •  K:-  ■  1=82 , 

. .  . . .  MiiUffiiii 


ILTsL  ..  ■==”~^ 


CNJ ( J ) 

Nj 

CON! 

ct 

COP (K  »N } 

COUNT 

CPOP 

K 

cpp 

P' 

COI Jt I *J) 

CRO 

Rc 

CT  OP 

tJ 

CTP 

T 

CTV J ( J > 

Tv, 

CUOP 

u; 

CVJ(J) 

vi 

CWI t I ) 

Wj 

CW J(J) 

% 

CXE J(Jt 

xej 

CZi < I > 

7  • 

*•1 

DELY 

DELYC 

DELYP 

^START 

DENN(  J> 

DEVF 

DF10I I  ) 

AF-/((L'TO 

DKF  Ml) 

B*. 

EECCPJI J) 

EECCCpj) 

EECHJ (J) 

EECChj) 

ECCNUU  <  J 5 

EECC/iJ/T) 

C-4 

NUMBER  OF  VIBRATIONAL  LEVELS 
41 »3228«CR0frCT0P 

COEFFICIENT  OF  X**<N-1)  IN  g,  C,  POLYNOMIAL 
LINE  COUNTER 
REFERENCE  PRESSURE 
PRESSURE 

DIMENSIONLESS  RATE  OF  PRODUCTION 
GAS  CONSTANT,  1.98647  <  CAL /MOLE -°t<) 

REFERENCE  TEMPERATURE 

s  , 

temperature 

VIBRATIONAL  TEMPERATURE 
REFERENCE  VELOCITY 
VIBRATIONAL  COUPLING  FACTOR 

DENOTES  WHICH  SPECIES  PRODUCTION  RELATION  V0  USE 

fyfl/Tvj  -  \/Tl 

NOT  USED,  AVAILABLE  FOR  FURTHER  SPECIES  DESCRIPTION 
DENOTES  WHICH  SPECIES  PRODUCTION  RELATION  TO  USE 
R4JNSE-KUTTA-R I CHARDSON  FULL  STEP  SlZE 
RUNGE-KUTTA-RICMAROSON  UPDATED  STEP  SIZE 
STARTING  VALUE  OF  STEP  SIZE  FOR  R-K-R  INTEGRATION 
NUF.3ER  DENSITY  CF  SPECIES j 
RESIDUE  IN  A  R-K-R  INTERVAL 

CHANGE  IN  STANDARD  FREE  ENERGY  FOR  REACTION  l 
FORWARD  OR  BACKWARD  RATE  CONSTANT  COEFFICIENT 
ELECTRONIC  EXCITATION  CONTRIBUTION  FOR  CCPJIJ) 
ELECTRONIC  EXCITATION  CONTRIBUTION  FOR  StfJfUi 
ELECTRONIC  EXCITATION  CONTRIBUTION  FOR  SVJOUT(J) 


EOF  2 

EpF  ;» 

i 

EPSJC J) 
EPSJ*.4<  j> 
ETA J ( J  > 
EXTRA J ( J  ) 
GAMMAO 
GJOP(J) 

I BETA  I ( I ) 
!BET!J( I « J) 
I  BOP 

ICAIJU  «  J) 
ICON 
IOELXC 
I  DUMP 
TEXT  ! 1 
IFAIL 
I  GO 
I  !  SF 
!  ISrPI 
!  ISG 
! NOSUM 
iOOP(K ) 
f  POT 
IRKIND 
I  RUN 
ISC 


INDICATOR  FOR  WRITING  END  OF  FILE  ON  TAPE  2 

INDICATOR  FOR  WRITING  END  OF  FILE  ON  TAPE  3 

€j  VIBRATIONAL  ENERGY 

fcjo,  VIBRATIONAL  ENERGY  AT  LOCAL  TRANSLATIONAL  TEMPERATURE 

rjj  NUMBER  OF  ATOMS  PER  MOLECULE 

TEST  j  TEST  PAR/ METER*  j=  !«•*** IS 

NOT  USED*  AVAILABLE  FOR  FURTHER  SPECIFICATION 

■x'  REFERENCE  SPECIES  CONCENTRATION 

J° 

?«•  P*  ~Z-  ?LJ 

frj 

INDICATOR  FOR  BOUNDARY  CONDITION 
Ay  VIBRATION-DISSOCIATION  COUPLING  INDICATOR 

NOT  USED*  AVAILABLE  FOR  FURTHER  SPECIFICATION 
S.S.C.  SUCCESSFUL  STEP  COUNTER 

DuMPjkd  INDICATOR  for  optional  OUTPUT  ON  TAPE  2 

TEST  it  (FIXED  POINT) 

MAIN  PROGRAM  FAILURE  INDICATOR 
STARTING  ADDRESS  OF  THIS  PROGRAM 
•C  DEFINE  VIBRATIONAL  NONCOUI LIBR  I UM  RANGE 

•f  +  I  ISF  4-  1 

g  DEFINE  VIBRATIONAL  NONEGU I L I 8~ I UM  RANGE 

E.E.  INDICATOR  FOR  ELECTRONIC  EXCITATION 

ORDER  OF  BOUNDARY  CONDITION  POLYNOMIAL 
INDICATOR  FOR  B.  C*  POLYNOMIAL  OR  TABLE 
R^RiW.  INDICATOR  FOR  THE  FOUR  LOCATIONS  IN  A  R-K-R  INTERVAL 
RUN  TOENTIFICATION 
C  NUM8ER  OF  ELEMENTS 


NMMMHM 


!SCP| 


c+! 


ISC  4  1 


s=>f- 

ISFP1 

IS(j 

ISGPl 

I  SR 

ISS 

I  START 

ISTOP 

IYST 

KFf INDt  J ) 


REND  IK ) 
REST (< } 

C-6 


•4  +  3+g  -  4- 


DEFINE  VIBRATIONAL  NONEQUfL I6RIUM  RANGE 
ISF  4  I 

DEFINE  VIBRATIONAL  NONEOUILIBR! UN  RANGE 
ISG  4  I 

NUMBER  OF  REACTIONS 
NUMBER  OF  SPECIES 
FINISHED  NUMBER  OF  SETS  OF  DATA 
NUMBER  OF  SETS  OF  DATA 
STARTING  STEP  INDICATOR 

INDICATOR  FOP  DIRECT! v <  CP  REACTION  RATE  CONSTANT 
INDEX  FOR  LOCATING  ELEMENTS  IN  B  AND  C  A  EYS 


Ml 

M+l 

INDEX  FOR 

LOCATING 

ELEMENTS 

IN 

B 

AND 

C 

ARREYS 

MSUMJfJ) 

mj 

NUMBER  OF 

ELECTRONIC  LEVELS 

•  MAXIMUM 

OF  0 

MX 

s+s-f 

INDEX  FOR 

LOCATING 

ELEMENTS 

IN 

B 

AND 

C 

ARREYS 

MX  1 

MX+  ' 

INDEX  FOR 

LOCATING 

ELEMENTS 

IN 

B 

AND 

C 

APprvjg 

MX2 

MX +  2 

INDEX  FOR 

locat ing 

ELEMENTS 

I.M 

B 

AND 

C 

ARREYS 

MX3 

MX  +  3 

INDEX  FOR 

LOCATING 

ELEMENTS 

IN 

9 

and 

c 

ARREYS 

MX4 

MX +4 

INDEX  FOR 

LOCATING 

ELEMENTS 

IN 

B 

AND 

c 

ARREYS 

MX5 

MX +  5 

INDEX  FOR 

LOCATING 

ELEMENTS 

IN 

B 

AND 

c 

ARREYS 

NC 

INDICATOR 

FOR  FULL 

OR  HALF 

INTEGRATION 

STEP  SIZE 

NOR 

NUMBER  OF 

B,  C«  POLYNOMIALS 

USED* 

MAXIMUM  OF  3 

NUf  1 1 > 

hi 

ij 

NU!  J(I  *U; 

LEFT-HAND 

SIDE  STOICHIOMETRIC  COEFFICIENTS 

HVlJPtUJy 

RIGHT-HAND  SIDE  STOICMIOMETRIC 

COEFFICIENTS 

ENDING  VALUE  OF  «  REGION  COVERED  BY  9*  C*  POLY,  -ft 
STARTING  VALUE  OF  A  REGION  COVERED  BY  B*  C«  POLY,  -ft 


RHOOF 


REFERENCE  JENS! TY 


ft' 


i 


R 


SAJCJI 

aj 

S3  J  ?  J  * 

bJ 

SGJL I J*L  J 

9ji 

SHJ(J) 

hj 

SHUOIJ) 

hj 

SHJOPI J) 

SKBI ( I ) 

Vi 

SKFf  (  I  ) 

V, 

SKFI  INU  > 

SPECIKI I *2) 

SPECJKI J.A) 

SUJOOTC J) 

fi/T 

TAUA^tU) 

S' 

TAUB  )ij) 

**S 

TAUCJIJJ 

TAUOU(J) 

TAUJP{J) 

xf- 

TEMP 

T 

THEVJ(J) 

9/j 

THEVJP ( J ) 

ei 

TP( J  00 ! 

tscale 

TYUOO! 

XA 

xe 

xe 

*j  -  bj  + 

CONSTANT  FOR  OSHICml  POTENTIAL 
DEGENERACY  OF  ELECTRONIC  LEVEL  Jt 
ENTHALPY  INCLUDING  MEAT  OF  FORMATION 
HEAT  OF  FORMATION 
HEAT  OF  FORMATION 
BACKWARD  REACTION  RATE  CONSTANT 
FORWARD  REACTION  RATE  CONSTANT 

FORWARD  REACTION  RATE  CONSTANT  AT  VIBRATIONAL  EGUIL 
REACTION  IDENTIFICATION 
SPECIES  IDENTIFICATION 
CHEMICAL  POTENTIAL  OF  SPECIES 

DESCRIBES  VIBRATIONAL  RELAXATION  TIME*  TAUJP(J) 
DESCRIBES  VIBRATIONAL  RELAXATION  TIME*  TAUJP(J) 
DESCRIBES  VIBRATIONAL  RELAXATION  TIME.  TAUJPIJ) 
DESCRIBES  VIBRATIONAL  RELAXATION  TIME*  TAUJPIJ) 
VIBRATIONAL  RELAXATION  TIME  OF  SPECIES 

temperature 

CHARACTERISTIC  VIBRATIONAL  TEMPERATURE 
CHARACTERISTIC  VIBRATIONAL  TEMPERATURE 
FUNCTION  VALUES  IN  THE  TABLE  OF  B*  C. 

CU0P**2*CMW0P/ < «  *  1 850 1 4E?*CR0*CT0P  » 
y  VALUES  IN  THE  TABLE  OF  S.  C. 

0.5#TSCALE 
6 , 022E23 /t CR0*C  TOP } 

CMWOP*CLP/ICUOP*RHOOP ) 

C-? 


X! ASIA!  XIA2(i>*U  ACCUMULAT ! ON  COUNTER  FOR  THE  R-K-R  INCREMENTAL  VALUES 

f> 

>nIA2(3>  *p 

XlA2  (4)  -  A 

XIA3CJJ  -gj  jsjr,.ys  ACCUMULATION  COUNTER  FOR  THE  R-K-R  INCREMENTAL  VALUfS 

XIA4(J)  £j  j*f+l,.~,g  ACCUMU*  ATfON  COUNTER  FOR  THE  R-K-R  INCREMENTAL  VALUES 

XLAMUfUI  VIBRATIONAL  RELAXATION  DISTANCE 

XLCT  inIT) 

S 

xNutm  vi  v\  i  v'i\ 

XNUI  J(  I  «  J)  2/-tj  LEFT-HAND  SIDE  STOICHIOMETRIC  COEFFICIENTS 

XNUIJPUtJJ  Vfj  RIGHT-HAND  SIDE  STOICHIOMETRIC  COEFFICIENTS 

VPREV  LAST  VALUE  OF  y  FOR  WHICH  RESULTS  WERE  PRINTED 

Y5T0P  y ST0p  YSTOPP/CLP 

YSTOPP  y^ajp  VALUE  OF  STREAMWISE  DISTANCE  THAT  WILL  TERMINATE  RUN 


Reads  in  and  writes  out  all  input  data  on  tape  #4. 
S4  Subroutine  WRITE  1 


V/ rites  headings  and  output  formats  on  tape  #4 
S5  Subroutine  WRITE  2 


Writes  results  on  tape  #4.  If  sense  switch  6  is  on  it  prints 
3 Cl  in  each  paragraph  on  tape  #4.  If  sense  switch  5  is  on  it 
prints  Qjj  on  tape  #3. 

S6  Subroutine  MATRIX 


Initializes  the  matrix  for  solving  simultaneous  linear  equations. 
S7  Subroutine  PER 

Computes  the  derivatives  of  every  variable  with  respect  to  y  . 
S3  Subroutine  SUB  1 

Writes  optional  outputs*  on  tape  #2 
$9  Subroutine  SUB  2 

th 

Computes  vibrational  relaxation  time  of  the  j  species  in 
the  l  "-*k  5  rsng®- 


Program  Listings 


Phis  section  presents  tht  complete  Port ran  IV  listings  of  tut 


mam  program,  as  well  as  all  the  subroutines  of  the  ftreamtubs 


program.  The  functions  of  all  the  subroutines  are  described  as 


follows: 


S2  Subroutine  START 


Performs  initialization  of  variables  and  definition  of  constants 


for  the  Runge-Kutta-Richardsen  method. 


Iter'- 


S3  Subroutine  READ  IN 


SI 0  Subroutine  3UB4 


mJ 

Computes  the  forward  rats  constant  py 
SI  1  Subroutine  SUB  5 

Writes  optional  outputs  on  taps  #2,  if  sense  switch  1  is  on. 

512  Subroutine  SUBS 

Computes  the  electron  excitation  contributions  to  the  thermo¬ 
chemical  properties  ilj  f  . 

513  Subroutine  TEST  1 

Performs  output  printing  step  size  test,  and  also  makes 
suitable  adjustments  on 

5 14  Subroutine  TEST  2 

f 

Performs  the  tests  for  negative  Y'  f  ~T  and  A 

51 5  Subroutine  TESTS. 

Performs  residue  test  for  the  Runge-Kutta-Richardson  method. 
If  the  test  fails,  it  halves  the  integration  step  size. 

516  Subroutine  OVERFLW 

Calls  subroutine  FPT  which  is  a  library  subroutine  to  skip  the 
floating-point  overflow  trap. 

Sll  Subroutine  LCOUNT 

Tape  #.'■  ejector 


518  Subroutine  STOP 

Performs  computing  termination  test  for  each  set  of  data  in 
this  run; 

519  Subroutine  SIMSO.L 

Solves  the  set  of  simultaneous  linear  equations  by  Grauss*  met7  >d* 

C-10 


ar 


■H*.  flj 


maa  w  * t  _  — 

o&v  *  aacuuu  xyxuxv  x, 

computes  derivatives  of  boundary  condition  fissctioa  wiib  respect 
tc  atreamwise  distance#  y  #  from  the  soundary  condition  table  by 
Foward  Gregory**Newton  formula# 

521  Function  YINT 

Interpolates  boundary  condition  function  from  the  boundary 
condition  table  by  second-order  Newton  divided-differences  method, 

522  Common  and  Dimension  Statements 

Used  with  the  main  program  and  all  tho  subprograms  except  sub 


routine  SIMSOL#  function  YINT  and  function  DERT, 


ISTAifT*© 

ASSIGN  3000  TO  ISO 

| ftftrtn  er*Sr 

3 SO I  CAUL  READ IN 
|YST*© 

3002  II5PP1*!SF*1 

3003  USG»?SG 
5004  CAUL  «?R|TE  I 

IffXTl  I “EXTRA J (Ut 
3020  1SCP1 *ISC+1 
3030  1SFP1 * ISF+1 
5040  ISGP1  » I  SG-f  1 
IISF  =  1SF 

5240  M*4+!SS+ISG-!SF 
3230  M1*M+1 
3260  MX*ISG-ISF4!SS 
5270  «Xl«MX+i 
5280  MX2*MX+2 
5290  MX3*«X+3 

5300  MX4*MX+4 
MX5«MX45 

5301  CALL  MATRIX 
5310  OELY»D£LYP/CLP 

YSTOPbYSTOPP/CLP 
CRO* 1*98647 

3330  TSCAlEb  <  CUOP*CUOP*CMWOP >/ ( CR0*CT0P*4 • 1 650 1 4E+7  > 

XC* ( CMWOP*CLP  > / <  CU0P*RH00P 1 
I0£LXC=0 

XB*6* 022E+23/ ( CR0*CTOP 1 

XA  *  *5*CUQP*CU0P  *  CMWOP/l 4*1 8501 4C+07  *CRO  *  CT0P1 
YPREV  *  -1*0 

5340  00  5370  J* 1*ISS«S 

5350  THEVJ( J?*THEVJP< JJ  /CTO P 
5360  SH JO ( J } =  SHJOP ( J >  /<CR0*CT0P1 

5370  SAJ(Jl*S8J(J}+*3*<5*0+  2*0* (ETAJ ( J 1-1 *0 1 1*AL0G(CTQP1 
5360  00  5400  !>:«!SC*1 

3390  DO  5400  J*1  « ISS* 1 

5400  8ETAI J  <  1  *  J  >  *XNUI  .IP  U  «  J  1-XNUI  Jit  •  J  J 

5410  00  5460  J» 1  *  1  SR* 1 

5420  XNUI ( I  1*0*0 

5430  ©ETA1 {f)*0*0 

5440  00  5460  J*1*ISS«1 

5450  XNUK  1  )*XNU1  ( !  >  +XNU1 J  ( I  *  J 1 

3460  SETA  I  ( I  1  *BETA  HI  H-BETAJ  J{  I,  J) 

7000  00  7050  J*i«lSS*i 

7010  DO  7020  L*1 *8*1 

7020  CE JLP ( J»L 1 -CEJLPX ( J«L  1/ <3*29820E-24*CT0P > 

7030  EECHJ1 J 1*0*0 

7040  EECNUJ(J1*0*0 

7050  EECCPJU  1*0*0 

7060  CONI* CRO *CT0P*4*  1 3228E401 

7070  DO  7120  l«!«ISR.l 

7080  00  7100  J*1«ISS«1 


50-1 

Nil 

501 

Sol 

SOI 

98! 

SO* 

set 

SOi 

SSI 

SOI 

SOS 

SOI 

SOI 

soi 

SOI 

SOI 

SOI 

ten 

soi 

SOI 

SOS 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOI 

SOS 

SOI 

SOI 

SOI 


23 

24  * 
23 

26 

a?  , 

25 

29 

30 

31 

32 

33 

34 

35 

36 


37 

39 


39 

40 

sj. 

42 

43 
C4 
43 

46 

47 
40 

49 

50 

51 


SOS  52  * 

SOI  S3 


C-12 


7100 
7110 
f  1  SO 

5501 

5010 

561 1 

5612 


5619 

5660 

5661 

5670 

5671 

5680 

5681 
5690 
5700 
5710 
571! 
5720 
5613 


5723 

5724 


5614 

5721 

5730 

5740 

5750 

5760 

5770 

5780 

5790 

5791 

5792 

5793 
5800 


5810 

5820 


ierri  jij*  ji-betaiju.js 

I8ETAII I )«8ETAI t I » 

NUf  1  1  >*XNU|  { I  » 

A21<41»CA0P 

IF  <  OELY  -  EXTRAJ<14>  >5611  ♦  56U  *1*1 9 

»RlTE<6,56ig>  I  RUN  «  A  l  T  51  IvDELY 

FORMAT <6X*4HRUN  * 46*9H  SK IPPC0.5X6HAT  Y 

ISTART*fSTART+l 

00  TO  5000 

00  5670  J«*l,3.l 

A! T! J )»At ! (J> 

Ai !C( JiaAl I (J) 

AgfCl JlaA2I (J) 

ASTiJi^ASt (J) 

A2!C(4)bA2! (4) 

A2TU>«A2!  C4> 

00  5690  J*I t I SS* 1 

A31C( J)»A3f (Ji 
A3T <  J)*  A3 1 ( J ) 

IFi ISG-tSPP* )  5613*5710 *5710 
DO  5720  J* I SFR1 *  ISO* 1 
A4ICC J>«A4I ( J) 

A4T< J)»A41 (J) 

NC«I 

IPUPOT.EQ*  0)  GO  TO  3614 

00  5723  J*1.N0R 

1F1AITU  WLT.TY(J)  >  GO  TO  5724 

CONTINUE 

IRAl.j 

tFdPAUEQ*  1>  IRAt-2 
IFURA1  «GE*NOR)  |RAl*NOR-l 
A2T(  1B0P)»YINT(7Y*TP*  |RA1  «AITU  >  > 
OELYC-OELY/FLOAT ( NC > 

DO  6240  1C«1 »NC 
DO  5740  J* 1*4.1 
YIm2( J>«0*9 
DO  5770  J«l*20«l 
XI A3( J>*0,0 
XIA4<J>»0*0 

DO  6130  l*l» 1*4*1 

IRK IND* 111 
IFAIL  »  0 
CALL  TEST2 

IF < IFAIL >  561 0*580© *5610 
CALL  DER 

IF CNC *NE* l  «GR«  lCaNE.l)  GO  TO  5820 

IFA1L*0 

CALL  TE5T1 

IF( IFAIL  >  5610*5820*5610 
IFCHl  «NE#  1>  GO  TO  5809 
IF1IC  «NE*  NC>  GO  TO  5809 
00  5801  *3 

41  1C(  U )*AI T ( J > 


*EI8*e*5X4«0Y  *£18*8} 


881 

ms 

SOI 

ms 

«i 

it 

m 

#4 

ms 

m 

SOI 

m 

sol 

SOI 

m 

8ft 

SOI 

8# 

SOI 

70 

$01 

71 

SOI 

78 

SOI 

73 

sot 

74 

SOI 

79 

SOI 

76 

SOI 

77 

SOI 

78 

soi 

79 

Sol 

80 

sot 

91 

SOI 

88 

SOI 

83 

SOI 

SOI 

89 

86 

SOI 

87 

101 

88 

SOI 

89 

SOS 

SOI 

90 

m 

so* 

SOI 

92 

93 

sot 

94 

SOI 

95 

SGI 

SOI 

96 

97 

SO! 

98 

SOI 

99 

SOI 

100 

SOI 

101 

SOI 

102 

SOI 

103 

SOI 

104 

sot 

105 

SOI 

106 

C  5,1 


3801  A2ICC jiasgfCJJ 
A2fCl4)'*A2T(4f 
DO  3802  J*i.ISj 

3802  A3IC’ J)»A3T« J) 

IP<  ( 130-ISffPl  )  *LT«  0  >  60  TO  5809 
DO  3803  J«ISFP1*ISG 
5803  A4 iC ( J  5 *  A4 T  ( J } 

5Q09  JJ=0 

DELYC*D£LV/FLOAT(NC  > 

5630  IFUISG-I  ISFP1  )  5930.5840.5840 
5840  00  5920  J*  JfSFPI  «  i  1SG.  I 
5850  1 

5860  TEMPI  «0ELYC*8<  JJ.MU 
5870  XtA4<  J)*X|A4CJJ+At©H(M  III  )#TEMPI 
5890  A4TC  Jl=A4iCI  JM-8ETA(  i 1  I >*T£MP1 
5920  CONTINUE 
S930  TEMP2»0*0 

00  6010  J*1 .ISS.l 

5940  JJ®!1SG-I ISF+J 
5930  TEMPI *DcLYC»B( JJ.Ml) 

5960  X|A3<  J)*:XIA3<  Jl+ACPHAt  I  1 1  )#TEMP1 
5980  A3T< J)«A3IC(J)*3ETA  (HI  )*T£MP1 
6010  T£MP2«  TEMPS  +  A3T?J> 

6020  DO  6100  J* 1.4*1 

6030  JJ=MX5-J 

6040  TEMP l »DEl_YC*S(  JJ.Ml  > 

6050  XIA2C  J)aXIA2(JH-ALPHA(  1  1 1  )#T2MP1 
6070  A2T< J?«A2!C( J)+8ETA< 1 1  I |*T£MP1 
6100  CONTINUE 

6110  AJTC1  1“A1  ICU  l+BETAUII  )*DELYC 

6120  A1T(31=1*0/TEMP2 

6121  A1T<2)»A2TI3>#A1T<3  M»TSCALE/m2T (2 ) 
6130  CONTINUE 

6140  IPmSG-nSPPl  >  6170.6150.6150 

6150  00  6160  J-1ISFP1.IISG.1 

6160  A4T(JJ»A4IC(J>+XlA4(J>/6*0 

6170  TEMP2*0,0 

6180  00  6200  J»1.1SS«1 

6190  A3TC  Jl*A31C<  J14-XIA3  (J  1/6*0 

6200  TEMP2-TEMP2+A3T<J) 

6210  DO  6220  J* 1.4.1 
6220  A2T(Jl«A2ICCJ)+XlA2(J>/6*0 
6230  A1T(3  5*U0/TEMP2 

AS  7iS )=A2T{3)#AlT(3J*rSCALE/A2T(2) 

6240  CONTI NUfe 

6241  GO  TO  J6242.62471 .  NC 

6242  NC“2 

6243  DO  6244  J*l,3 
AJTI ( JJ*AIT(JJ 
A2T1 ( J}«A2T(J) 

A 1  T  C  J )  *  A 1  !  f  J  > 

6244  A2TC J?*A21 (J) 

A2TJ  <4  J*A2T(A 1 


A2TC4|aA2!44i 

life 

Do  6245  J*-  i  *  j  SS 

SOI 

itr 

A3T1  ( JJ*A3T(U» 

sot 

ISM 

&S45 

A3T;\J;«A3I (Ji 

sot 

163 

IFttl  ISG-I  !SFPn*LT*0)  <30  TO  5614 

sot 

tm 

00  6246  J*  1  f SFP1 *  l  I SG 

SOS 

163 

A4T! { J)344T( J) 

SOI 

6246 

A4T( J)»A4I iJi 

SOI 

tot 

GO  TO  5614 

sot 

16# 

6247 

1F( < 1 ISG-I 15FP1 1 ,LT.O)  GO  TO  6250 

Sot 

160 

ifail«o 

SOI 

170 

CALL  TE5T3<A*T*A4T1« JISFP1.I1SG  *4 > 

SOI 

i7i 

JF( IFAIL#NE«0J  GO  TO  5610 

SOI 

172 

00  6248  J* 1 ISFP1 • I ISG 

SOI 

173 

6248 

A4T ( J ) *A4T  ( J)+(A4T(J)-A4T1 ( J) )/l5# 

SOI 

174 

6250 

TEMP2*0#0 

SOI 

175 

f  FA I  L*0 

SOI 

176 

CALL  TE5T34 A2"sA3Tl « 1 « !SS»3 ) 

SO! 

177 

IFUFA!L.NE*0>  GO  TO  5610 

so; 

1 78 

DO  6249  J»1 ♦ ISS 

SC! 

179 

A3T(J)=A3T  ( J)+(A3t< J1-A3T1 < J) >/15* 

SOI 

180 

TEMP2»TEKP2+A3T(J> 

sot 

131 

6249 

continue 

SOI 

132 

I FA 1L*0 

SOI 

183 

CALL  TEST3IA2T.A2T1 .1 .4^2) 

sot 

134 

IFUFAIL  *NE.  0)  GO  TO  5610 

SOI 

189 

00  625!  J«! «4 

SOI 

isa 

6251 

A2T( J)*A2T(J)+CA2T( J)-«2T|  (  J)  )/|5* 

SOI 

187 

AlT(3l»l#0/TEKP2 

SOI 

183 

A! T{2 ) *A2T { 3  t*A 1 T (3 )*TSCALE/A2T  <  2  > 

SO! 

t  AC 

6260 

GO  TO  5613 

SOI 

190 

END 

SO! 

191 

SUBROUTINE  START 

i  mM 

CALL  AVRFLK 

SOS 

* 

1000 

00  1040  J*l»3«! 

30aJ 

$  m 

1010 

sot 

€ 

1020 

A1T<J1*0*0 

sot 

3 

1030 

A2f <31*0*0 

&ot 

0 

1040 

A2T1 J)»0*0 

SOS 

m 

f 

1041 

A21 (4  1*0*0 

aoi 

9 

1042 

A2T (4 1 *0*0 

SOE 

0 

10S0 

00  1090  J*1 *20* 1 

SOS 

so 

1060 

A3i CJ1*0«0 

SOS 

ft 

1070 

A3T  C  J ) *0*0 

SOS 

ifi 

1080 

A4!< 31*0*0 

SOS 

IS 

1090 

A4T<31*0»0 

SOS 

14 

1190 

EOF2*0*0 

SOS 

ss 

1  110 

£OF3*0.0 

S02 

to 

1  120 

ALPHA <1 }®I*0 

SOS 

I? 

1130 

ALPHA (2)*2*0 

SOS 

16 

1140 

ALPHA (31*2*0 

SOS 

S9 

1 1  SO 

ALPHA  <  4  J  *  1 ,0 

SOS 

go 

1  160 

BtfTACll*  *3 

SOS 

k'l 

1  170 

BETA! 2l**5 

S02 

t« 

1160 

BETA! 31*1*0 

S02 

S3 

1  190 

BETA! 4 1*1 *0 

SOS 

SO 

1200 

return 

SOS 

29 

END 

SOS 

26 

* 


C-I6 


SUBROUTINE  READ IN 

01  N£NS ION  ISGJL(B).  NALP1J120).  1PP14) 

INTEGER  00G1  s  0002*  QQ34t  OnMt  G*SG?.  PfP! «  PPPSt  l§3 

ALPHAMERIC  WORDS  IN  INTEGER  FORM  SG3 

0002  *  COMP, 

QQQ2  a— 17508747779  f$3 

Q6G4  *  NO  563 

QQ04  1 7997957478 
0005  »  YES 

0005  .-17997969234  S.03 

GOG  *>  .FORW  SS3 

0006  *  24270826544  S®3 

Q0Q7  .REV  S03 

0007  a-J H030C91 31 2  ®03 

IPPU  )»  VEL#  S®3 

SPPi  I  )•*- 1807471 1792  ^3 

IPP(25»  RMO  S"3 

I PP ( 2 1 ■- I 7S74 1 85264  503 

|PP(3 )* PRESS,  S®3 

IPP<3 ) .—8209771675  S03 

I PP ( 4 ) «  AREA  S03 

IPP(4).  17475916912  S03 

112;  REAO  1 130. IftUN* ISS. ISR. ISF* 1SG» f Sc* IBOP* INDSUM. 1 DUMP.DELYP. YSTGPP  S03 
1130  FORMAT  <  A6. 8 13«  2E1 0,4  J  S03 

1140  I7EA0  1  ISOtCPOP.RMOOP.CTOP.GAMMAO.CACP.CUOP.CMWOP.CLP  S03 

1150  FORMAT  (8E9.5)  503 

1160  READ  M  70. ( EXTRA J ( J ) «  J*  f  *  1 S *  1  >  s03 

1170  FORMAT  ( 5E 14,7)  33 

1180  DO  '  230  J*1«1SS«1  r'03 

1190  REAO  !200*ETAJ( J) •  SB J ( J ) » THEV JP ( J ) ♦ SM JOP( J ) . CNJ ( J ) • CXEU( J  )  S03 

1200  FORMAT <©E1 2.6 »  s®3 

1210  REAO  1 220  «  TAUA J( J ) »TAUBU( J ) • TAUCU ( J  > TAUOJI J ) t 1 SPEC JK i J ) oK* I ,4 1503 
1220  FORMAT <4EI 2»6» 4A6  )  "  S03 

1230  READ  1 240 « MSUM J ( J ) • ( SG JL (J.L)*L*1*8»1  ) »  1 CE J  <( J.L) .L»l .6. 1  )  S03 

1240  FORMAT! 12. 8F2.0/8E9. 5)  S03 

1270  DO  1300  )*1 .1SR.1  S03 

1280  READ  !290fCWI (I).CZI(l )*CDI < I) ♦ <XNU! JPC I . J > . J*1 «20 1  *  S03 

1  (XNUlJ;i*U)«Jal *20) • ICAIUt I , J>«J*i»20>  S03 

1290  FORMAT  ( 3F2»0 . 2GF 1 , 0 « 20FI ,  0 » 20F 1,0)  S03 

13CC  READ  13J0.KFI !ND( 1) »AKF! ( I) .BKF! U ' *CKFI ( ! l.DKFI < 1  1  « ISPECtK U «K> ♦  S03 

1 K* 1 *2 . 1  )  s03 

1310  FORMAT ( I4«4E!4»-.?A6'  503 

1320  00  1330  U=l»ISC»;  S03 

1330  READ  1 340 o  I ALP! J i I  * J > ♦ ! * 1 * 1 S<» 1  )  !  03 

1340  FORMAT  <20F2,0)  SO 3 

1 SFP1 * I SF+ 1  S03 

1350  IFUSG-ISFPl)  1370*  1360*  1360  S03 

1360  READ  1170,  tA4MJ>»  J=1SFP1«1SG)  507 

1370  READ  1170.  (A3I(J>*  JM.1SC)  303 

READ  1170.  A!  I  U  )  .A!!  (2  1.A21  (3)  .A2112  ,*A2!  ( I)  *A|  I  <3*  SO’’ 

;  JPOT«0  -  POLY.  1P0T»1  -  TABLE  S03 

REAO  1378. NOP. IPOT  S03 

1370  FORMAT  1213 J  SO 3 


5  ip§p|: 

m 

103 

m 

$03 

S7 

S#3 

is 

IF  1 1  POT  *N£  •  01  TO  1385 

t~79  do  1382  IRal.NQR  -  5iT^ 

READ  13S0,  JC«P<  IRI.RE3TI  !R).RSND<  IRjt  S©3 

1380  F0RMATU3.2E14.7)  §03 

13B1  lOOPi  !R)«IOOP< 1R)+I  fl09 

I OP* I  OOP { f  R  J  SOS 

1382  READ  1170.  <COP<!R.J}»  J-1.I0P)  $83 

GO  TO  2000  J03 

1385  READ  1386.  'TY< ]R) ,TP< IR) *  lR*i .NOR)  S03 

1386  FORMAT < 2E *  4*7)  ©03 

2000  WRITE  ( 4<;  2001  >  ©03 

2001  FORMAT (6H1  RUN. 5X. 1HS. 5X«  I  KR.5X .  1 HF . 5X. 1 HS. 5X.  1 HC .2X14B3QUN0ARV  C$03 

1 OND. . 3X.2HEE . 8X . 5H0ELYP , 1 OX . 6HYS  fOPR . 6X . 4K0UMP )  $03 

2008  IFUNDSUW)  201 1  .2009,20 1 1  $03 

2009  PPP2*0Q04  $03 

2010  GO  TO  2012  3©3 

2011  PPP2*  OCK35  $03 

2012  WRITE  54.2013) IRON. ISS. ISR, ISF. ISG, ISC, IPP< 130P1 .PPP2.DELYP.YST0PPS03 

!  .IDUMP  503 

2013  FORMAT < 1 X. A6* 15*4 16.6X, A6.4X, A6.2X,  1  P2E1 5.7.4X.I1 )  $03 

2014  WRITE  (4*2015)  $03 

2015  FORMAT ( 1H0.4X.4HCP0P. 1 1  >  JHRHOOP, 1 1 X.4HCT0P. 1  OX.  $03 

1  4NCJOP.IOX.5HCMWOP.1 1X.3HCLP)  S03 

2016  WRITE  (4.201 7 )CPOP. RHOOP.CTOP.  CUQP.CMWOP.CLP  $03 

2017  FORMAT (1P8E1 5.7)  $03 

2016  WRITE  (4.2019) (EXTRACT  J  >  .J*  1 .  5,  1)  $03 

2019  FORMAT ( S  5H0  CONTROLS  AND, 1 P5E1 5.7 )  $03 

2020  WRITE  <4.2021 ) (EXTRAJ ( J ) . J®6. 10*1 )  S03 

2021  FORMAT <15H  CONSTANTS  .1P5E15.7)  $03 

WRITE  (4,30005 (EVTRAJ(J),  Jail, 15,1)  $03 

3000  FORMAT  < 15X.1P5E15.7)  S03 

2022  WRITE  (4.2023)  $03 

2023  FORMAT ( 1H0 , 1 IX, 1 H J, 8X.4HETA J« 1 1 X.3HSB J* 11 X.6HTMEV JP.9X. 3HSHJ0P* 1 1 XS03 

I . 3HCN J , i 2X .4HCX0 J )  S03 

COUNT* 11.0  $03 

2024  DO  2028  J*l,J2S.l  S03 

WRITE  (4*2025) J.ETAJ(J) .S8J(J).TMEVJP(J) *SMJOP(J) ,CNJ(J)»CXEJ{J)  S03 

2025  FORMAT (1 IX. 12, 2X, 1P6E15.7)  S03 

2026  COUNT »CGUNT+ 1.0  $03 

2027  CALL  LCOUNT  $03 

2026  CONTINUE  $03 

2029  WRITE  (4,20.10)  SQ3 

£030  FORMA  T  (1  HO  «  1 1 X  ,  1 H  J  .  7X  ♦  5r, ,  AU  A  J  ,  1  OX  »  5HTAU3  J  ,  1  OX  «  5HT  AUC  J  •  1  OX  ,  5HT  AUD  J  ,  SG3 

113X«7H*<=*EC!ES)  S03 

2031  COUNT *COUNT+2.0  $03 

2032  CALL  LCOUNT  $03 

2033  DO  2038  J-l.ISS.l  $03 

2034  WRITE  <4,2035 )J, TAUAJIJ ) ,TAU8J< J) ,TAUCJ( J) .TAUOJ ( J ) *  (SPECJKS03 

1 (J.K) .K>1«4.1 )  303 

2035  FORMAT < 1 1 X, I2.2X. 1P4E15.7.4A6 )  $03 

2036 . COUNT  «COUNT+ 1 . 0  $03 

2037  CALL  LCOUNT  $03 

2038  CONTINUE  $03 


C-1S 


ocejip 

7 


WI  Zl 

16 


at  nui  jp 

CAI J  ! 


1 1 
16) 


2039  WRITE  (4*2040 

gC.40  FORMAT < I 20H0  J  M  SSJ!  234367  SCEJ1P  s 

j  4  5  6  7 

2041  COUNT “C0UNT+2.0 

2042  CALL  LCOUNT 

2043  DO  2030  J»1*1SS*1 

2044  00  2043  L*i«8*l 

2043  IS6JL(L)aSSJL( J«L> 

2046  WRITE  (4*2047)  J*MSUMJ(J)«  <  ISGUL(L)«L«1  **L*I  )  * 

j  (CEJLPXt J*L)*L*I«S«I > 

2047  FORMAT (1 3* I2*2X*8I3«F5* 1 ♦ IP7E12.S) 

2048  COUNT "C0UNT4-1  «0 

2049  CALL  LCOUNT 

2050  CON“INUE 

2059  WRITE  (4*2060) 

2060  FORMAT (I06H0  I  WI  Zt  DI  NUIJR  I  6  11  1© 

IIJ  136  11  16  CAI J  1  6  11  16) 

2061  COUNT «COUNT+2#0 

2062  00  2066  I*1*1SR*1 

2063  00  2066  J»!«!SS«1 

2064  NUIJP(I*U)*XNU!JP(1 «J) 

2065  NU1J(I*U>  ■XNUIJ(I«J) 

2066  1CAIU ( I « J) »CA! J( I • J ) 

2067  DO  2073  I«)*1SR*1 

2068  NWI*CWltI) 

2069  NZI-CZI ( I ) 

2070  NO I “CO I ( I ) 

2071  WRITE  (4*2072) I «NWI «NZI *NOI • (NUI UPC I *J) »U*1 *20  *11* 

1 l » J) « J*1 *20  • 1  ) « i ICAi J( i * J ) * U»i «  20*1) 

2072  FORMAT (3X* I2«2X« I2*2X« I2*2X« I2.9X.S1 1 « IX *51 1 • 1X*5I 1 • 1X*3I 
1  IX*5I  1  *  1 X  «  5 1 1  *  1 X  *  5 1 1  *  6X  ♦  3 1 1  «  l  X  •  S 1 1  '■  1 X  «  5 1 1  *  1 X  •  5 1 1  > 

2073  COUNTsCOUNT+1 *0 

2074  CALL  LCOUNT 
2073  CONTINUE 

2076  WRITE  (4.2077) 

2077  FORMAT (1 SHOO I RECT SON  I  .6X.4HAKF1 « 1 1 X.4H8KFI • 1 1X.4HOCFI 
1 4H0XF I ♦ 1  OX . 8HRCACT 1 ON ) 

2078  COUNT  « COUNT +2  * 0 

2079  CALL  LCOUNT 

2080  DO  2089  I*I»ISR*1 

2081  IF  ( KF 1 1  NO  (  D)  2084  «  2082  *2084 

2082  PPPl»OGQ6 

2083  60  TO  2085 

2084  PPP1 ■QQQ7 

20 es  WRITE  (4*2086)PPPI« I*  AKFf ( I ) *BKFI ( I ) «CKFI ( I J «DKF! < I )  * 
IK(I«K)«K»I *2*1 ) 

2086  FORMAT (3X*AS*2X* 12* 2X* l P4EI 5*7* 3X» 2A6 ) 

2087  COUNT  »C0UNT+1 *0 
2008  CALL  LCOUNT 

2089  CONTINUE 

2090  WRITE  (4*2091 ) 

2091  FORMAT (lOHOALPJX  J«1 *7X» IH5*8X»2H10*8X*2H15*8X*2H20) 

2092  COUNT *C0UN?+2*0 


I-  3*3 

)503 

SOS 

$03 

503 


«6X*4HAKFI « 1 IX*  4H8KP I » 1 1 X  *  4HCXF 1 


803 

SOS 

S03 

303 

SOS 

NUS03 

SOS 

SOS 

S03 

$03 

$03 

$03 

§03 

$03 

$03 

SOS 

$03 

(NUIJ(S03 

$03 

1.9X.SI1.S03 

$03 

$03 

§03 

503 

$03 

«!1X*  SOS 
$03 
$03 
$03 
$03 
$03 
$03 
$03 
$03 
(SPSC1S03 
$03 
$03 
$03 
$03 
$03 
$03 
$03 
503 


110  ■ 

111 

m 

ns 

m 

ns 

no 

it? 

ns 

no 

tto 

121 

122 

123 

124 
2SS 
126 
127 
t2S 
120 

130 

131 

132 

133 

134 
139 
136 
1ST 
139 

139 

140 

141 

142 

143 

144 
149 

146 

147 
146 

149 
190 
tSI 
192 
153 
194 

155 

156 
197 

150 
I  199 


019 


■  vi  i 


£393  CALL  LCOUNT 

f094  00  2i0f  J«iOSC.l 

2095  DO  2096  I*l«|SS*l 

2096  NALP! J( ! )*ALPIJ( I* J) 

2097  WRITE  ( 4  »  2093 1  J*  C  NALP |J(1)* 1*1*1 SS •  I  ) 

2098  FQRMATC4X, I?*2X*20|2) 

2099  COUNT *COUNT+ 1 *0 

2100  CALL  LCOUNT 

2101  CONTINUE 

A2 1 (4 ) »CA0P/CLP*»2 
2200  WRITE  (4*2210) 

2210  FORMAT (/19H01NI T I AL  CONDITIONS/I  HO *6X IHY* 14X1 HA*  14X1HT* 14X1HP* I3X 
1  3HRH0* 13X1 HU* 13X2HMW) 

WRITE  (4*2017)  A1 I ( 1 )*A2I (4 )♦ A1 I (2 >*A£ I (3)* A2I (2) • A2I (1 »«A1 1 (3) 
COUNT  *COUNT+6,0 
CALL  LCOUNT 

WRITE  <4*2220 )  ( A3I  ( J ) * J*1 *  I SS > 

2220  FORMAT ( 1K0 *58X2HGJ/ < 1 P8E15.7 >  > 

COUNT  *COUNT+5,0 
CALL  LCOUNT 

1FUSG-ISFP1  )  2250*2230*2230 
2230  WRITE  (4*2240)  < A4 I < J ) * J* ISFP1 • I SG ) 

2240  FORMAT ( 1H0 .57X4HEPS J/( 1 P8E1 5«7) ) 

COUNT *COUNT+5,0 

CALL  LCOUNT 

2250  IF( IPOT.NE.O)  GO  TO  2260 
WRITE  (4*2241)  IPP(IBOP) 

2241  FORMAT ( 1H0 . 15X5HRANGE *32X8HC0PU  OF  »A6) 

COUNT  sCOUNT +2 • 0 

DO  2251  IR»1«N0R 
IOP*IOOP( IR) 

WRITE  (4.2242)  REST ( IR) *P£ND ( IR ) * (COP ( !R* J) ♦ J*1 * t OP) 

2242  FORMAT  ( 1PE15«>7  <3M  TOE  1 6. 7.  I  H*  *6E  1  5o 7 ./ ( 33X*  1  P6E15.7  )  ) 

EL  IN* I00P( lR)/6+l 

COUNT  *  COUNT+EL I N 
CALL  LCOUNT 

2251  CONTINUE 
GO  TO  2262 

2260  WRITE (4*2261 )  |PP(fBOP)«  (TY( IR ) *TP( IR) *  !R«J.NOR) 

2261  FORMAT (1 HO «40X*26HT ABLE  OF  BOUNDARY  COND*  -  «A6«6H  VS,  Y#/ 

1  ( 1 X2E1 6,8* I X2E16,8« 1 X2E1 6,8*1 X2EI6«8) ) 

EL  IN* NOR/4+3 
COUNT *C0UNT+EL IN 
CALL  LCOUNT 

2262  RETURN 
END 


m3  %m 


2250 


2251 


2261 


3-5$  tm 
S03  164 
SOS  34$ 

sss  i  m 
ses  .6? 

503  86* 
£03  J69 
593  170 
$03  1TI 
303  172 
SOS  ITS 
SOS  174 
SOS  ITS 
S93  ITS 
503  IT9 
503  178 
$03  179 
503  too 
503  101 
$03  182 
SO 3  183 
S03  184 
503  163 
$03  186 
$03  187 
503  188 
503  189 
S03  190 
$03  I 91 
SOS  192 
$03  193 
$03  194 
S03  195 
503  196 
$03  197 
$03  193 
S03  199 
$03  200 
$03  201 
SOS  302 
S03  203 
S03  204 
$03  1105 


C-20 


suoromtine  imitEi 
DIN©*$!0*«  6(20) 

INTEGER  21  *  23*  24*  25*  26*  Z?«  28*  Z9*  8U«C»  XI*  X3«  X3*  X«  * 
XS*  XA*  XT*  XS*  X9*  XI 6*  XII*  XI2*  XI3*  X14*  X15.  IH6* 
XI7*  XI8*  XI9*  X20*6 

zi*  §Am& 

21 *—1767035137 » 

23*  EPS 
23*- 17997647026 
24*  TVJ 
24«»I79979707£n 
25*  CHI 
ZS*- l 7997837849 
26®  G 

26*-! 7997956! 64 
27*  P 

27*- 179979581 63 
28*  RHG 
Z8*-1 7997927974 
29*  VEL 
29*-! 7997976752 
BLNK* 

6LNK*-1 799795SJ92 

XI  TO  X20  *  1  TO  20*  RESPECTIVELY 

XI *-17209429000 

X^-- 172262 06256 

X3®- l 7242963472 

X4*-l 7259760668 

X5*- 17276537904 

X6*-l 729331 31 20 

X7*— 1731 009233C 

X0*-I 7326669552 

X9«- 1 73^/3646766 

XI 0*1066524464 

XI 1*1 103301680 

XI 2* 1120078696 

X!  3*  U  36656ns 

XI 4* 1153632326 

XI 5*1 1 7041 0546 

XI 6* 1187187760 

XI 7* 1203964976 

XI 8* 12207421 92 

XI 9* I 23751940a 

X20*2I602662«S 

G(l)  *  XI 

GI2)  *  X2 

6(31  *  X3 

6(41  «  X4 

6(5!  *  X5 

@(6>  *  X6 

G<7>  *  X7 

G*©>  *  XS 

6(9)  *  X9 


$04  It 

m 

$66 

S04  4 

*a>* 

-9 

& 

§06 

l 

S@4 

$04 

f 

$06 

to 

$04 

it 

$04 

t& 

$04 

IS 

$04 

16 

$04 

15 

$04 

16 

$04 

17 

S04 

2$ 

$04 

19 

$04 

20 

$04 

21 

$04 

22 

$04  23 

$04 

24 

$04 

£3 

$04 

26 

$04 

27 

$04 

26 

S04 

29 

304 

30 

S04 

31 

$04 

32 

$04 

33 

$04 

34 

S04 

35 

$04 

36 

S04 

37 

S04 

35 

$04 

39 

$04 

40 

504 

41 

S04 

42 

$04 

43 

$04 

44 

$04 

45 

$04 

46 

S04 

47 

$04 

49 

S04 

49 

$04 

56 

$04 

51 

$04 

52 

$04 

S3 

ones*  xio 
0(11)*  xu 
eua»«  X12 
6(13)*  XI 3 
6U4h  X14 

ansi-  Xi2 

<S(16>*  xi6 
OUT)*  X17 

sub)*  a ia 

6(19>«  X19 
G(2C)«  X20 
ISCP1  *  1SC4-I 
COUNTaO,0 

100  PRINT  101 « IRUN 

1 01  FORM* T (4H1 RUN* A6.2X  » 1 HY  *  1 3X.4H0ELY  *9X*  26HT 
IE ) 

1 02  WRITE  C4.lo3)IRUN 

103  FORMAT C 1HI »7X*4HRUN  *A6) 

104  WRITE  (4*105) 

105  FORMAT  ( 1 HO • 50X * J ©HFORMAT  FOR  RESULTS) 

1090  WRITE  (4*11001 

1 1  OOOFORMAT  (8H  Y*  1 4X*  1  HT«  I4X*  I  HA*  14X*  ll)P*  1 3X*3HRH0*  I3X*  1HU*  1 3X* 

I  2HMW) 

1110  WRITE  (4*ll20)(2l*I*I«l*ISS*l) 

1120  FORMAT  (8(  3X*A6*  12*4X  >/8( 3X *  A6«  I2«4X )/4  ( 3X «A6«  1 2*4X )  ) 

WR1TE(4*10) 

10  FORMAT (40X*42HNUMBER  DENSITY  PER  C*  C*  p0R  EACH  SPECIES  ) 

CALL  SSWTCH (6*K000FX ) 

(SO  TO(2O0O*20O2)*KOOOFX 

2000  WRITE  (4*2001 > (Z5*I *I«1 *ISR*1  ) 

2001  FORMAT (8 (2X«A6* I2»5X)/S(2X*A6* 1 2 *5X >/8 (2X* A6« !2» 5X )/8(2X* A6« I2*3X )S04 

1/8(2X*A6*I2*5X>  >  304 

2002  IFU  1S6-I ISFP!  )  106*1150*1150  304 

1  ISO  WRITE  (4*1 160) (Z3* ! *  I«I ISFP1 *  1  ISO* 1 )  304 

1160  F0RMAT(8(2X»A6*12*5X)/3(fcX*A6«r2*5X)/4(2X*A6* 12*3X)  )  S04 

1170  WRITE  (4.1 180) (Z4* I •  !» 1 1SFP1 o 1 1 SG. 1 )  304 

1180  FORMAT  (8(2X.A6* I2.5X >/8(2X. A6. l2*3X )/4(2X. A6* I2.5X)  )  504 

106  WRITE  (4.107)  304 

107  FORMAT (//////SIX* 16HFCRMAT  FOR  DUMPS)  304 

108  WRITE  (4.109)  S04 

109  FORMAT! I HO *6X* 1HY* 1 4X» 1 HT* 1 4X*2MMW* 1 3X* 1  HU* 1 3X.3HRH0* 1 3X* IMP* 12X»  S04 

13HRUN*4X*5KRK1N0>  s®4 

UO  WRITE  (4*111  )(Zl*J.I»l*  JSS*  1  )  304 

111  FORMAT (8 (3X*A8* I2.4X )/8 (3X» A6« I 2.4X) /4 ( 3X. A6* I 2*4X ) )  S04 

|F(  1 1  S6**l  1 SFP1  )  130*  112*112  304 

112  WRITE  (4*S13)(Z2*I*  1*I1SFP1 . 1ISG* 1 )  SO* 

113  F0RMAT(8(2X*A6.I2.3X>/8(2X.A6.l2*3X)/4(2X*Ae. I2.5X))  S04 

130  DO  114  J*1 « ISS* 1  304 

114  WRITE  (4*115) (U*Kal *8*1 )  304 

115  F0f®<AT<10H  SUJOOT* I2.8X.4HCCPJ* I2*9X*4HEPSU*I2*9X*3NSHJ* I2*9X*6S04 

lHEPSJIN*!2*ax*4HCTVJ« I 2 * 8X » 3HXLAMU . 12 .9X.3HCVJ* 1 2 )  504 

116  00  117  I*1*1SR*1  304 

117  WRITE  (4* 1 18) ( I *K*1 *7.1 )  304 


504  63 

S04  64 

504  65 

504  66 

504  67 

RATIO  OR  TYPE  OF  FAILURSC4  66 

504  69 

504  70 

504  71 

504  72 

504  73 

504  74 

504  75 

504  76 

504  77 

504  76 

304  79 

504  SO 
504  81 


U8  v0as$AT(3X«7HD€LT  FI  *l£«tXi#HOd*!*12f  18S»3H0KI  f_0T 

1  7X»4hSkF3  *I2*9X»4HSKBI*I2*9X*4HCH?  !  *  125  #|| 

119  150  120  1*1  e  ISO*  1  SSflS 

*20  WRITE  (4*121  >  (Z6*  1*0*  0* * 5CP1 • I  S3* 1)  ,Sg  v-|J[ 

121  FORMAT  (80X.A6*  12*  12*2X1  )  f?4  ®|l. 

122  WRITS  (4*  5 £3) 

123  FORMAT! 2 20H  (PARTI  A(_S  OF  VARIABLES  LISTED  BELOW  WRT  ¥  1#| 

1  6/LINE  -rjm  lit 

124  IF(  I ISG— I ISFPi  )127*  125*  123 

125  WRITE  (4* 126) (Z3*6< I ) • 1*1 ISFPI *  1 BS6* l > * (21 *3( I )• t*l* ISS  «  i  I  «ZT$(?Uiti@#  II© 

1K*ZB*BLNK*Z9*6LNK  |*T 

126  FORMAT  (8(3X*A6«A6) >  S2* 

50  TO  126  *  1^ 

127  WRITE  <4* 126? (Zl *5( 1 1 « 1*1 • !SS*1 ) «Z7*BLNK*Z3*8LNK*29«ILRK  SO#  ISO 

128  WRITE  (4*103) IRON  *?* 

RETURN  164  t&M 


sutsmitra^  mite  2 
DIMENSION  dshn  <  20  i 

1007  WRITE  (4*  1009  JA1 S  <  I  ;«A.i  ic-!)v*;;f  =  }«A11<3) 

i  COS  FORMAT  C///1P7E15*©) 

1009  write  (4«ioio) cAat(J)*J«*t«tss*n 

1010  FORMAT  UP8£1S«6) 

DO  10  J*t*  tSS 

10  D£NM<  J)e6*0a30eS3#A31  ( Ji#Ag$  ( 2  >  fc^OOP/CSWOP 
WmTE(4.10I0>  (0£t#J(J)«J»!«!*£) 

CALL  SSWYCH ( 6  *  KQOOFX ) 

SO  70(3000*3001 >*KCOOFX 

3000  WRITE  U*!010)(Cmi(I>*!*«i«ISR*i> 

3001  IFUISG-II5FP1  )  3004*3002*3003 

3002  WRITE  (4*1010)  (A4I  (U),  Jtaf  fSFPl  *  !  ISO*  !  ) 

3003  WRITE  (4*1010) (CTVJ ( J) *  1 ISFP1 • I  ISO* 1  ) 

3004  CALL  SSWTCH(5*K0G0PX) 

SO  TO(3OOS*3O1O)*KO0OPX 

3005  WRITE  < 3*3006 )A1 i(l)« I RUN 

3006  FORMAT (///6X*1P£15*6.9X«A6) 

3007  CO  3008  !I*i«lSR*! 

3008  WRITE  (3*1010) (CO! J ( 1 1 • J) • J*I SCP1 *  I SS* 1 ) 

3009  £QF3*1«Q 

3010  COUWT*COUNT+I*0 

3011  IF ( COUNT-EXTRA J( 1 >)  3015*3012*3012 

3012  WRITE  (4*3013 ) I RUN 

3013  FGRMATdgHl  RUN  ,45) 

3014  C0UNT*0„0 

3015  CONTI NUE 

3018  CALL  STOP 
YPREV  «*  A1  I  U  ) 

3019  RETURN 
END 


SOS  9 
SOS  IQ 
S0s?  1 1 
SOS  12 
SOS  13 
$8S  £4 

SOS  IS 
SOS  ;© 


SOS  17 
SOS  I© 
SOS  10 
Sss  20 
SOS  21 
305  28 

SOS  23 
SOS  £4 
SOS  25 
SOS  2© 

**•**•?  e*  » 

SC5  £8 
SOS  29 
SOS  30 
SOS  31 
SOS  32 


subroutine 
teoo  oo  1  cos 

10©J  DO  1002  LL®1  #5^1*1 

1002  €(L*U.)»0«0 

1003  ««»! 

1004  ir<n$G~i isFPi  Hoostioessioos 
ioos  oo  loo?  j*i isppi  « use? i 

1006  C<MS*UMtf)*l  .0 

1007  HN«WN+1 

1008  00  1012  Jmi ♦ 1SC« 1 

1009  DO  1011 

10|0  JJJim  )}  ISG-I  itF*JJ 

1011  C(K«»JJJl)a 

1012  mm«MM+1 

1053  DO  1015  J»ISCPl *1SS#J 
1014  C<f«*WWM»«l ,o 
|0l5  MMaf?M  +  l 

JP»MX3-1B0P 
CIMXl »JP)*>1 oO 
1017  C<MX2»KX2»»I«0 
101©  C ( NX3  «  WX3 ) « J  ,  0 
1019  RETURN 


sueaourtNB  ©es 

c  kodsl  for  stpeak  tum 

too©  call  sum 

S0©|  IF  UWSUSH003*  1003*1002 
;oog  CALL  SUSA 
toes  Tg®»*Air<2j 
IOIO  KLCTpAL€K?<T£'MP) 

1020  CTP»CTQP*Tg&P 

1030  CPP  a  A2T(31*CU0P4CUOP*»H0O» 

10*0  JFilSFl  1200ei200«:050 

1050  00  1 190  J«1 •ISFol 

106©  IF'CETAJC J1-I.01  1  12Q*  1070#  1120 

1 070  SUJO0T  <  J } =«SA J { J } -2  *5  *XLCT+SH JO « J ) /TEMP-EECNUJ < J 1 
1060  CCP  J  C  J  5-fEECCPJ « J  J 

1000  £PSJ(J>*0#O 

1 1 00  SH  J  <  J )  s2*5«TEMP*SHJ0  t  J  > 4EECH J  i  J I 

II  !0  GO  TO  H90 

1420  AA*EXP(STHSVJ* JJ/TEMP4 

1130  AAitt#S*(5e042#0  #  (ETAJ 1 J )- 1 ©0 >  ) 

1140  AA2*  £TAJ<  JJ  — 1 *0 

1 1  SO  SU JOOT  t J l m-SAJ 5  J 1 -A A 1 #XLCT+ AA2* AL06 1 \ .0- 1 • 0/A A 1 4SH  JO  C  J I /TE«p- 

1  EECNUJ(J» 

1160  gPSJf  JS*TH£VJ(J>/{AA-l,01 

1170  CCPJ< Jl*AA!+AA2*EPSJt JJ**2*AA/TEMP»*2  4£ECCPJ<Ji 
1160  SHJ<J)sAAl*TEMP+AA2«EPSJi JJ4SMJ0CJ14EECHJI J> 

1150  continue 

1200  IF  USS-iSGPl)  1300 »1210«I21O 
1210  00  1350  J*ISGP1 *!5S*i 

1220  IFCETAJI J1-I.05  t 280* J230e 1 280 

1 2  30  SU JOOT  C  J ) » -SA J i J 1 -2  «5*XLCT+SH JO  <  J ) /TEMP-EECNU J  *  J  J 
1240  CCPJl J>«2*S+EECCPJ(J) 

1250  EP5J<J»«0«»0 

1 260  SH  J  (  J  >  «g.5*TENP*SH JO  ( J 1 -frESCK J  <  J 1 

1270  GO  TO  1350 

1260  AA*EXP ( THE VJ  t  J  > /TCNP ) 

1200  AA1«  *5  *  C5*0  42*0  *  (ETAJ ( J)-l .0 )  1 
1300  AA2*ETAJ»:Ji-l»0 

1 31  0  SUJOOT  C  J)*=~SAJ  * J 1 -AAl #XLCT+AA2*AL0S(  1  *0-1  * O/AA  1+ShJO  C J  )/T £WP- 
1  EECWUJU) 

1320  EPSJI J)»TH£VJ(J1  /  (AA-1.0) 

1330  CCPJC J!*AA14AA2*EPSJ( J>*&2*AA/TE«P#*2  +EECCPJ(J> 

1340  SHJ« J)*AA1«TEMP  4AA2^EPSJ ij)  +SHJO ? J 1+EECH J( J J 
1350  CONTINUE 

1360  IF  C 1 SG— I SFP1 >  1560* I  361  * J 36  < 

1361  CALL  SU32 

1370  00  1540  J®  I SpPl * l SG* 1 

1360  B8*EXP*TH£VJ( Jl^TEMP) 

B3l*  ETAJ«J)-1*0 
3®2«  ,5#(5*042*e#B@! 1 
1300  EPSJiNC J>*THFTv'J?JJ/<B8-l  ©0) 

140©  CTVJC  JlssTHEVJf  Jl/ALOGl  (THSVJC  J14A4TI  J)  )/A«t«  J)J 
1410  CCPJCJ1«®82  +EECCPJ ( J } 

1420  SUJOOT  C  Jls-SAJ?  J1-8B24XLCT+B®! 443.0 6<  1  .0-1*0/60  l+SH  JO  f  J  >/T£NP~ 
C-26 


3B*h  S 
SOYA  f 
SOYA  3 


$074  6 
YOtA  7 
SOYA  S 
SOYA  f 
SOYA  10 
S07A  11 
S07A  1:2 
SOYA  13 
SOYA  14 
S07A  IS 
SOYA  16 
SOYA  IY 
SOYA  IS 
SOYA  10 
SOYA  gO 
SOYA  21 
SOYA  22 
SOYA  23 
SOYA  24 
S07A  25 
507A  26 
SOYA  27 
S07A  28 
SC 7 A  20 
S07A  30 
SOYA  31 
SOYA  32 
SOYA  33 
SOYA  34 
S07A  35 
SOYA  36 

Soya  37 
SO 7 A  3S 
S07A  39 
S07A  40 
SOYA  41 
SOYA  42 
SOYA  43 
SOYA  44 
S07A  45 
397A  46 
SOYA  47 
SOYA  40 
SOYA  49 
507A  5© 
S07A  51 
SOYA  52 
SOYA  53 


£tm  m 


LOWER  BOUND  OF  CM! i 


<737  CONTINUE  S$?AI6f 

TEMPI®  C5M  <tl  1*TEMP1#(RH00P*A2T<2)  I?I  I  +  GDI  (If  I4TEMP1*  S07*HSf* 

I  (RMO0P*ACT(2>  >**2/(AiT(3>*C«t(0Pl  807*309 

1740  TE«P2»  C£I ( ! I >*TENP2*(RHO0P*A2T(2» >»*2  /CMW0P**2  S07AtJO  „ 

1738  DO  1739  J*!SCP1 • tSS. l  S07AI 1 ! 

1739  CO  I  J<  1 1*  J)  *XC*BETAI J(  1 1  *v»)*SKFt  (  1 1  )*CHt  I  Cl  I  )*  (TEMPI  +TgMP2*A?T<  Jt  jS07A$t* 

1790  CONTINUE  &Q7AH3 

1800  00  1820  L»!cM*l  w07AU4i 

1010  DO  1020  UL«t*MI#I  S07AII5 

1820  B(L«LL>®  CCL.LL)  S©7AU* 

BCMXl *Ml )*0*0  S07A1I7 

I F  ( I POT  4  NE  #01  GC  TO  1870  50 7 A 1 f  0 

IRA»I  507AH9 

GO  TO  (1820*1027*18 25)#  NOP  S07A120 

1826  IF«A1 T(1 >*GT#REND(2 > )  IRA* I RA41  S0TAI21 

1827  IF(A1T(1  ) * C«T •  REND ( 1  )  1  IRA»IRA41  S07AI22 

1828  MN«I00P( IRA1-1  S07A123 

IF(MN>  1830*1030*1823  S07A184 

1823  DO  1829  J«1*MN  50TA12S 

WN1»IOOP( IRAl-J+l  5O7AI20 


FMM*MNl— 1 

1829  8S«X1 *Ml J»8(MX1 *M1 )  *Al~(l )+FMN*C0P( fRA«MNt ) 

GO  TO  96 

1870  DO  1871  J«I*NOR 

'F(A1T(  1  )*i.T.TV(J))  GO  TO  1872 

1871  CONTINUE 
1072  JRA»J 

IF(IRA*E0*1)  IRA»2 

1F( 1RA*GE*N0R)  IRA*NOR-l 

B»MX1 *M1 >»OERT(TY*TP* IRA*A1T(1 )) 

06  CMLL  5SWTCH ( 4  *  KSW ) 

GO  TO  (97*I83n>*KSW 

97  KR 1 TE  ( 6  « 98  5  A1TC 1 > • IRA«B(MX1 *M1 } *DELY 

9©  F0RNATC4H  Y  *E16.8*6X*5HIRA  *14«0X4MDF  »EI  6.8.6X.4HC  f  *-£16.0) 
1630  0CMX2.WX4J*  A2T(2)*A2T < 1 ) 

1840  NN»i 

1350  IF  (IS6-ISFPIJ  1690*1860*1850 

1860  DO  1 880  J» I ISFPl *  I I3G* I 

TEMPI *0*0 

T£MP2»O*0 

DO  1861  I  I “1  *  1 SR. 1 

TgMP3*(CAI J(I 1 *J»*CQI J( 1 1* J) )/  (  A3T( J}#A2T(2 )*A2T( 1 )*CHI I ( 1 1 ) ) 
TEMPI "TEMPI +TE MP3 

1861  TEWP2-TEMP2+TEMC3* ( 1 *0-CHI 1 ( 1 1 >  > 

TEMP3«THEVJ(J)*U  #0/CVVJ( U>-1  .O/TEMP) 

Y£MP4« (EPS JIN ( J1—A4T I J ) )/XLAMJ( J) 

TEMPS®  ( THE V J (J)/( EXP ( TEMP3 ) - 1  * 0  1  *>(^NJ(  J  >*THEVJ(  J )  )  / 

1  (EXP(CNJ( J)*TEMP3)-U0)  -  A4T(./>  )  *TEMP1 
TEMP2  *  -(*S»(CNJ( J)-1*0)*THEVU( J)  *4T(J>)«  TENP2 
S(NN*«U  «  TEMP4+TEMP5+TEWP2 
3880  NN»NN4l 

189o  DO  1940  J«ISCP1 * 1SS*1 
1900  JJ«1 ISG-1 5SF+J 


S07A12Y 

SOTAIRfi 

S07A129 

S07.*l  30  " 

S07AI31 

S07A132 

S074133  n 

S07A134 

SC7A13S 

S07A136 

S07A137 

S07A138 

S07A139 

S07A140 

S07A14I 

SG7A142 

S07AI43 

S07A144 

S07A145 

SO  Ml 46 

S07A147 

S07A148 

S07A149 

S07A1S0 

S07A1S1 

S07A1S2 

S07A133 

S07A154 

S07A155  ' 

S07AIS6 

S07AI87 

S07AS56 

S07A139  ' 


C-28 


1910 

5920 

1930 

1940 


Tff  MP1  9  0 

00  1930  !  *  I • f  SW  « l 

f E«rja  TEMP!  4-  caSJU»3! 

S(3J*M1  i*  TEN"  /  <A2Tf21*  AST  C  t  1  5 


*.  Ml  mm  »  _.  j\  A 

;jf  Tr:  m*-  I  ^ 


i960 
l  ■<-' 70 
1900 
*990 
3000 
201  0 
2020 
2030 

2049 

2050 
2060 
2070 
2080 
2090 
2100 
2101 


2105 

*  !  0 
'0 
30 


OO  1 970  3=1 « IS5* 1 

tfMPI  *  TEi^3j  +  A3T<  J)*CCPJi  3) 

TEi»Pls  TEMP»*  TSCALE 

p{?dX4  ^XSle  TEMP!  *Ai  T 1 3  5  /&2T121 

0{MV.4tMx3i®  -  ?  TEMPI  *A!  T  <3 1*A2T  13)  5  /A2Tf2?4*S 

aCMX4,MX4>*  TSCALE  *  A2TU3 

LLL*0 

IF ( ! 1 SG“I I SFP1 12070* 2040» 2040 

00  2060  3* 1 !SFP! *  1 ISG* 1 

LLL»LLL*S 

B(MX4,LUL1  *«3TI31*{£v43C J1-! *0) 

KISS*  ISS 

DO  2  5  '-'0  J=!«I33;jl 

KXKaJ  JSG-I  kSF+3 

rj  i  MX4  « K*tX 1  a  SM3<  31- (TEMPI  *4  !  T  <  3 }  **2<*A2T  <  3 !  >  /A2T12) 
B 1 MX3  »MX4  >  *A£T  121/A2T  <1  1 
B(H  ;3*MX1 1*A2TC2?/A2T<4 i 
CALL  SSWTCH<3*N5W3) 

fF(NSH3*EQ*l 1  PRINT  2105»AITC 1 1 »DSLY 
PORMAT <7H0AT  V  *Eia«B«  SX4H0Y  «E18«8> 

CALL  S1MS01.  C  B*M?45) 

CALL  5US5 

RETURN 

END 


sotAieo 

S07A16I 

S87AI6S 

5L7Ai63 

S07*1S4 

S07AIG3 

S0YAi6& 

507*167 
S07A168 
S07MS9 
507A 1 7© 

S07A 1 71 

SGTM7S 

SO 7 A 1 73 

S07A 5  74 

S07A l T5 

S07A176 

S07AI77 

507*178 

507A i 79 

SO7A10O 

S07AS8! 

S07A182 

S07A ! ©3 

$07*364 

507* 1 85 

S07A188 

* ?7A187 

5U7A1&8 


|pSjp 
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6002 

6006 


SUBROUTINE  ICR 

PREFERENTIAL  MOOEL  FOR  STREAM  TU®£ 

I F i YPREV  >6000*1000*1000 

6000  IFnSS-ISFPl?  I000«600l  *600  i 

6001  WRITE  (4*6007) 

6007  FORMAT (5X* 1KN*  ?X* ! HU* 1 1 X*2HWE *  1  I  X*  4HWEXS *  1 0X«4HWEY£  *  1 0Xs4HW£2ES 
00  60C4  J*  ISFPi « ISG* 1 

READ  6002*NTEMP< J).UTEMP(J)*WE( J)*WEXECJ>*W£YE<J>.WeZElJl 

>002  FORMAT  (12*  SE 14*7)  ^ 

WRITE  (4*6006W^JTEMP(J>*UTEMP(  J>*WE(J>  *WEX£(  J)*SCYE(U>«  W£2E<J1 

>006  “ORMAT ( !6* 1P5E14*7 ) 

NTEMP ( J ) *NTEMP ( J  >  + 1 
TEMP2*0*0 
temp3*g*o 

KKP  *  NTEMP (J) 

oo  6003  kkm.xkp.i 

XKK  *  KX 

EV-P( KK. *  J ) * • XXX ■  •5  >  *  ( WE  ( J )  +  ( XKK- *3  >  * < “WEXE ( J  >  +  (XKk  - #3  >  *  <  WEYE < J 1 
1  +WEZ£ ( J  >  *  ( XKK- *  5 )  )  > ) 

IF (KK— 1 )  6009*6008*6009 

5008  TEMP20-E7JPU  «  J? 

5009  EV JP( KK ♦ J )  *  EVJP(KK, J1-TEMP20 
TEMP4«EXP(EVJP(KK* J)/UTFMPt  J> ) 

1£MP2»TEMP2+EVJP(KK* J)*TEMP4 

6003  T£MP3sTEMP3+TEMP4 
«JBAR ( J  >- ( 1 .9Q7322  *TEM P2  >  /  < TEMP3*CT0P*CR0  > 

6004  QJMUJ( J)*T£«Pt 
♦  000  CALL  SUB I 

1001  IF (INOSUM) 1003* 1003*1002 

1002  CALL  SU86 

1003  TEMP* A 1 T (2 ) 

!  G!  0  XLCT*ALOG(TEMP> 

1020  CTP*C TOP* TEMP 

1030  CPP  a  A2  r  (3)*CUO**4CUOP*RHOOP 
1040  IFCISF)  1200*1200*1050 
1050  DO  2190  J* 1*1 SF . 1 
1060  IF (ET  AJ(J)-1«0)  1120*1070*1120 

1 070  SUJOOT ( J ) *-SA J ( J ) -2  «  5*XLCT+SH JO ( J ) /TEMP-EECNU J ( J  ) 

1080  CCP J ( J  >  *2 • 5+EECCP J ( J ) 

1090  EPSJC  J ) *0«  0 

l 100  SHJ(J)*2*5*TEMP+SHJ0( J ) +EECHJ ( J ) 

1150  GO  TO  1190 
1120  AA*EXP(THEVJ« J)/TEMP) 

1130  AAl*«5*(5.0+2e0  *  (E7AJ ( J )- 1 « 0 )  > 

1140  AA2*  ETAJ(  J  »-l  *0  ^ 

1  130  SUJC07  ( J)*-SAJ(  J>-AA1*XLCT4-AA2*AL0G<  1  *0- 1  •  0/&A  >  +  St .  JO  ( J  ) /T,_MP- 

1  EECNUJ(J) 

1160  FPSJ( J)*THEVJ( J)/« AA-1 *0) 

U70  CCRJI J)»AA1+AA2*EPSJ(J)**2*AA/TEMP**2  +EECCPJ<J> 

1  1 80  SH J ( J ) *AA 1 *TEMP+ AA2*EPS J  <  J ) +SHJO <  J  14EECHJ ( J ) 

I  190  CONTINUE 

1200  IF  (ISS-ISGPl)  1360. 1210* 12'iO 
1210  00  1350  Jj«  I  S6PI  *  ISS*I 


6008 

6009 


6003 

6004 
♦  000 
1001 
1002 
1003 
!0!0 
1020 
1030 
1040 
1050 
1060 
1070 
1080 
1090 
1  100 
1150 
1120 
1 1  30 
1 1  40 
1130 


SOT®  s 

SO*m  3 

§6781  # 

sots  s 
sots  6 

S07B  t 

sots  s 

SOTS  9 
sots  10 
sots  11 

$070  12 
SO70  13 
SO70  14 
SO70  15 
SO70  16 
SOTS  !  t 
SO"©  1® 
S07B  19 
SO 7©  20 
SOtB  PH 

sot©  1:2 

S07B  23 
SOtB  2* 
SOTS  25 
SOtB  26 
SO 78  27 
SOtB  2© 
SOTB  29 

5078  30 
S07B  31 

5079  32 
S076  33 
S078  34 
sore  35 
5078  36 
SO 78  37 
S078  30 
S07B  39 
S078  40 
SOtB  41 
SOtB  42 
SOtB  43 
SOT©  44 
SO70  43 
S078  46 
SOtB  47 
S07B  40 
S07B  49 
SOtB  50 
S078  31 
SOTB  S2 
SOT©  33 
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ItaO  !P<£TAJU>-1#0!  l£SO« 1230* 123® 

1830  SUJOOTt J)«-SAJ<J1-£,5SXLCT+SHJ0{ JJ/T£WP~geCNUJ< J> 

1540  CCPJt J1®2»5+£ECCPJ1 J} 

S250  EPS JtJi*Q0O 

I860  SHJ<J)»g*s*TEMP+SHJO( J>*££CHJ(J> 

1270  GO  TO  1350 

1200  AA»EXP1THEVJ{J)/7EMP> 

1290  AA1 »  *5  4  15*0  +2*0  #  1ETAJ1 j)-j ,0 )  ) 

1300  AA2«STAJ{Jl-l,0 

1310  SUJOOT  ( J )  a— SA  J  f  J  > -A A 1  *XLCT+ AA2*AL0G  ( 1  *  0- 1 . 0/  A  A I 4S  W  J0  ( J  5  /TEMP- 
1  CCCNUJ(J) 

1320  EPSJt  J)*TH£VJU)  /  {AA-i*OS 

1330  CvP  M J)»AA1+AA2*£PSJ(3S*»2*AA/TE«P*«2  <fEECGPJ<J> 

1340  SHJ(J).AA!*TE*P  +AA2*EPSJ<J>  +SHJ0 1 J 1+EECH J( J > 

1350  CONTINUE 

1360  :F  (ISG-1SFP1)  1560*1361.1361 

1361  CALL  SU82 

1370  DO  1540  Jm  | SEP] * lSG* 1 

I3©9  Ba«gXP(TM£VJ(J )/TE«P) 
eai«  ETAJ(ji-i*o 
882*  *5#<5*0+2#0»BBl > 

1390  EPSJf N(J)*THEVJ( J)/ (BO-1 *0 ) 

1 400  CTVJl J 1 «THEV J 1 J ) /A LOG ( < TM£VJ 1 J  >  + A«T  1 J  >  >/  A4 T  <  J  >  1 
1410  CCPJ(J)*B82  +EECCPJ 1 J } 

1 420  SUJOOT  ( J )  — SA  J  (J)  -BB2*XLCT+B8 1  *ALGG  11,0-1.  O/BB  14-SH  JO  1 J » /TEMP- 
1  EECNUJ1J1 

1430  CUJ1 J 1  *THEVJ 1 J }* 1 1 .O/CTVJl J)-l ,0/TEMP ) 

1440  SMJ1  J)o0S2*TEMP+BB1#A4T{J)+«5HJO1J)+EECJJ1J) 

1460  XI. A" J (J)m  TAUJP1  Jj#£uOP*A2T(  1  )/CLP 

TFJP1 J)»l*0/(1 *0/(CTVJ( J)*CT0PJ~1,0/CTP-1 .0/UTEWP1 J) i 

?EMP2*0.0 

?EMP3*0.0 

TE^4«0.0 

TEMP5-0.0 

KXP*NTEMP{J> 

DO  6005  KKal.KKP.i 

T£WP6*EXP  t -EV J® i KK • J » /TPJP (J) ) 

TEMP7*£XP{-EVJP(KK. J)  /  CTP  ) 

TEMPSbEXP  1  ~EV  JP 1 XX «  J }  /  (CTVJ1  J)*CT0P)  ) 

TEMP2aTEMP24-EV.}P  1 KK  •  J  )  *T£MP6 
TEMP3oTEMP3+  TEMP6 
TEMP4«  TEMP4  +  TEMP? 

6003  TE«P5»TEM9g  4-  TEMPS 

E@APJ{  jj«(  j  e9«*?3£2»TEMP2;  /  1 TEM  :->3*C90*C70P ) 

1470  1F1CWJ1  J1M520*1471  *1320 

1471  CVJ(J)=1S0 

1472  60  TO  1540 

1520  CVJl J}* (T£MP4*TEMP3 }/ 1 TEMPgs  QJMUJ1JJ} 

1540  CONTINUE 

1560  CALL  5tm 

1561  DO  1790  1 1* 1 v 1 5P» { 

1570  OF 501 1 1  1*6.0 

1 5t Q  DO  1 390  Jb 1 » ! SS » 1 


#07#  WB 


20  M  m 

S3.fi  £9 
SOTS  60 
SQ7S  61 
SOTS 
S078  63 
Sd?a  64 
SOTS  60 
SO 78  M 
SOTS  6? 
S©7©  6S 
SOTS  69 
SOTS  70 
S078  71 
SOT©  72 
SO  70  73 
SOT©  74 
8076  75 
S0TB  ?6 
SOT®  77 
S07B  76 
SOTS  79 
S07B  SO 
S07«  gi; 

5075  S2 
SOT®  @3 

SOTS  04 
SO70  85 
S078  66 
SOTS  87 
SOTS  80 
SOTS  89 
SOTS  90 
SOTS  9! 
SOTS  92 

5076  92 
SOT©  94 
S07S  93 
S07B  96 
SOTS  97 
307©  9® 
S073  99 
SOTS 100 
5076*01 
S07SIC2 
5079 103 
S078I04 
SOTS 105 

so  mi o& 


C«*31 


1590 
s  fir-o 

2010 

163© 

1640 

1630 

1660 

1670 

1680 

1690 

1700 

1701 

1702 

1703 
1710 
17!  1 

1719 

*720 


1 


3  723 
1726 


1730 

1731 

1732 

1733 

1734 


1733 

1736 

1737 

1 

1740 

1736 

1739 

1790 

1300 

ISSO 

1820 


Sri?!!!!!2rIouI,+seTA,J(j**j>#stwooT{jj 

"'i  * '  *  *  '“cKf'i-uf  io(in; 

riLLl  ft,*5KP!  C  1 1  %/i  £TEMP*eONl  >** IBCTAI ( I f  H 

f  is.  w*i  *i(iO 

JFUT©-J5FP1  5  1 670 *  1 650 *  1 630 
£30  1  660  J*fSFPJ  •  JS6*  1 

TgfcP!«TC*Pi#cVJI3>**!CAIJ(  jj*  Ji 
SKF|( | I»*SKFI 1N( H )  4T2WPJ 
SKBKfh  •  SKF1  N(|t)  /  CKKIH 

XXTEMP  w  T2HP1 
TEMPI  «U0 

00  1711  J«1«!SS.I 

IF  (  I  BETIJ(  l!«35)17D2*|  703*  1 702 

IF  (A3TUn  1710*1719*1710 

IF  CXNUIJI l i ,j)  51702*171 1*1702 

CONTINUE^1* 

GOTO 1720 
TEMPI «0*0 

CM  1 1  (  !!)*■(  ,0—f  (PHOOP  »A2T«2n*#isgTAI  {|t)5#TEMPl  / 

ipIcmmi7,°P?r  °T°’  exnAJ(s>  '  L0*sa  bound  of 

Ir f CM F IT C ! ( |  #£Q*  0*0)  GO  TO  1723 
IFtABSCCH!  m  n  ULT.EXTRAJISJ  >  CHfl  (  m-0»0 
GO  TO  1726 

^r,Um-LT‘  E^TPAJ<45>  CMIMIII.0.0 

TEMP2«0,0 

IF  (CZI (II >11732,1732,1730 
DO  1731  sl«!  *  ISS*  1 

IF  «CWI < II 151733*1733*1734 
IF (CD I (III  51740* 1740* 1734 
TEMPI »1 

DO  1737  J* 1  * { SS« 1 

IF  (A3T (J ) 5 l 736* i 732* | 736 

IF  <NUIJ(  n»J)H73$«n'37*1736 

TE«Pl«TEMPi  #  (A3T( J5/CMW0P)**NUIJ< n *3) 
continue  ,1W’ 

Tempi*.  CWI  (  m*TEMP14(9HO0P#A2T(25  5#«NU1  <n> 

( RKG0P&A2T (2  5  5  *#2/ ( A 1 T ( 3 ) *C MWOP ) 

TENP2*  C2I ( I ! >«T£MP24(RHOOP*A2TC2) J*#2  /CM WOP 
00  1739  J.ISCP1* JSS*I 

CONtS^J>  “XC#SETAI  J<  * 1  *  J,*SKFI  (  1 1  >*CH*  I  H  I  )*  (TEMPI 

00  1820  L«1 «M* 1 

00  1820  LL»I *M1 * j 

B (L*LL 5 «  C (L»LL 5 
8(MX1 ,MI }«0«0 

!P(IPQT*NE*0)  go  TO  1870 
IPAal 

GO  TO  (1028*1827*1826)* 


NOR 


C»32 


tOTffj * | 

limits 
S079U9 
56781 14 
S071IIS 
£078116 
30781 |7 
S07©f 18 
SO TBS  19 
SO 78 120 

507812s 
£078 122 
SO  78!  23 
3078124 
8078126 

(CKf  ( I  *  5*)(XTEMPS07S126 
3078 127 

CHH  SOTS i 28 

5078129 
3078 130 
SOTS 13J 
S07B132 
SOTS 133 
S078134 
S078 1 3S 
SOTS I 36 
S078137 
Sf/7SI38 
5078 139 
SOTS I 40 
SOTS 1 4  S 
307® 1 }g 
5070143 
SOTS 144 
S07814S 
S07B146 
SO70S47 
507814*1 
SO 76 149 
S07BI50 

♦TEMP2»A37(3) 55078*51 
5078132 
S078I33 
S07B154 
SOTS 8 65 
S078565 
SOT® 15? 
SOTS 158 
5078*59 


+  CO! ( 1 1 )*TEMPI* 
*2 


1829  a 


1870 

1B7| 

1872 


#AS  T < 1  )+pMN»CCPURA*NNl  ) 


60  TO  1 872 


182*  !F(A17Cn*GT*WENO(2H  J9A«IRA*! 

1327  jrCAlTCl  }ftGT*Rt.’®(l  »  >  fRAwJRA+1 

1828  MN»!OOP<|RAl-! 

|F(MN>  1830*1830* 1323 

1823  DO  1829  J»1*MN 

MN 1 » I  OOP  U  RA  >  -  J+ 1 
FMN»KNI~1 

1829  B(«Xt .Ml )«B<MXI *M! 1  *4 J T <1 )+FMN»eCP< IRA *MNf ) 

GO  TO  96 

1870  DO  1871  J*1«N0R 

!F(A1T<1 >*LT#TV<J> }  60  TO  1872 
1B7J  CONTINUE 
1872  IRA* J 

IF ( IRA*EQ» 1 1  1RA*2 
IFURA«GE*NOR)  !«A*NOR-I 
8  <  MX1  .Ml  )*DERT  CTY *  TP* JRA»AlT(  I  ) > 

96  CALL  SSWTCH<4*KSW) 

GO  TO  <97* 1 6301 *KSW 

97  WRITE  < 6*98  J  A1 T< 1 » « tRA.EHMXl *M1  J «D£LY 

98  FORMAT  C4H  Y  -2I6»8«6X«3HIRA  *14*eX4HDF  *£i 6*8*6X*4H0Y  *216.8) 

1830  B(MX2.MX4)»  A2T<2**A2T<1> 

1840  NN»! 

1850  IF  USG-lSPPl)  1890*1860*1860 
1860  DO  1880  J«J 1SFP1 • 1ISG* 1 
TEMPI *0»0 
TEMP2*0,0 

DO  1861  11*1*1  SR* 1 

Te«P3« <CA1 J< 1 ! *J)*C01J<!I*J))/  <  A3T ( J)*AgT< 2)«A2T< 1 J*CHf i < 1 1 
TEMPI *T£MP1+TEWP3 

1S61  TEMP2«TEMP2+TEMP3*(  1 ,0-CHf  I  <  1  11) 

XTEMP 1  *  <  EPS J! N <  J ) -A4T  <  J )) /XLAMJ <  J ) 

XTEMP2*(EBARJ( J)-A4T<J) >  *  TEMPI 
TEMP2*  —  tGJSAR (JI-A4T  < J ) )*TEMP2 
B(NN«M1 )  *  XTEMP1  +  XTEMP2  +  TEMP2 
1880  NN*NN+1 

1890  DO  1940  J*rsCP! * ISS* l 

1900  JJ*USG-I  15F+J 

1910  TEMPI *0*0 

1920  DO  1930  1*1 • 1  SR* 1 

1930  TEMPI*  TEMPI  +  CQUUtJ) 

1940  8CJJ.M1  )*  TEMPI/  <A2T<21*  A2T11M 

1950  TEMPI *0,0 

1960  DO  1970  J-J*1SS*1 

1970  TEMPI*  TEMPI  +  A3T < J )*CCPJ< J) 

1980  TEMPI*  TEMPI*  t SCALE 

1990  B<MX4«MX2>*  TEMP1*A1T<3)  /A2T<2) 

2000  8<MX4*MX31*  ~<*~EMPI*A1T  C3)*A2T(3>  )  /A2T'2>#*2 
2010  0  <MX4  *MX4 ) *  TSCALE  *  A2TJ1) 

2020  LLL»0 

2030  IF<  HSG-I  ISFP1  12070' 2040*2040 
2040  DO  2060  J*  V 1 SFPi *  I ISG* 1 
2050  LLL*LLL*1 

2080  B <MM4 *LiL )  »A3T C J)* <ET*  J< J >-l *0 ) 


A3T(J)*A2T<2)«a2T<1  i*CHil  (ID) 


507818! 
S87#l8t 
S878263 
SO 78 164 
£078! 65 
3070 I 64 
5070167 

sots 183 

SOTS 169 
SC7SI70 
S07BI71 
SOTS 172 
S07B173 
S07BI74 
S078IT5 
SOTS 1 7* 
S07B177 
S0?8 I 78 
sots 179 
SOtE »  30 
S07B 1 81 

sots  leg 
sots 183 

SOTS 1 84 

S07S 165 

S078I86 

3078 1  St 

S07B1B8 

3075189 

S07© 190 

S07S 191 

SOTS 192 

S07B193 

£076194 

S07BI95 

SOTS 196 

S078 1 97 

S07B196 

S078 199 

S07B200 

S07B201 

S078202 

£078203 

S076204 

S07B205 

S07S206 

S07S207 

S078208 

S07B209 

SO 78210 

S07B21 1 

3078212 


<1070  |SS 

2080  DO  SlOO  J»l«iSS«! 

foaa  seicsta  1 1  SG“  f  f  sf+ j 

£100  Bf«X4*KKtO*  SHU f J)~ (TEMPI *A* T <3 |4*2*Af T <3> 
2101  mMX3*MX4>«A£TC2)/A2T(U 
BtMX3,MXl )«A27(2)/A2T<4> 

CALL  S£3TCH(3«NSW3* 

!F<NSW3*£G#n  PS  I  NT  21 0S»  A I  T  U  )  *DELY 
21  OS  WOPNATITHOAT  V  *E18.0*  5X4MOV  *£18*81 

2110  CALL  S1MS0L  <  B«PT»4 5) 

2120  CALL  SUSS 
2130  RETURN 
0 


/AST (21 


wgfmww 

Mmmm 


S#7»2f7 

sotisi® 

S87S2I® 

SOtStgO 


S07B22S 
$079223 
SOT® 224 


SUBROUTINE  SUS5 

1000  CALL  SSWTCHl 1 .KCOOFX) 

SO  TO <100 1 .1008) *KO0Of« 

1001  irtlRS<lMD-IOU*tP)  1002*1002®  1003 

1002  WRITE  <2*|003)A1T<!  ? «  A1  T<2)  *  M7  <  3)  *  AST  <  1  i  »A2f  S  2  )  ♦  AST  i  3 )  *  AgT  <4  > 
1  »1RUN. IRKIND 

1003  F0RNAT</////lF'7EiS,7*3X.A6,  16) 

F0P2* 1*0 

1005  WRITE  <2.1006)<A3T< J).  J*>1«ISS«1) 

1006  FORMAT < 1 R0E1 5*71 
IFUISG-IISFRl  )  1008®  1007*  1007 

1007  WRITE  <2«100fe) <A4T<U)«JbIJSFRI. I ISS.l) 

1008  RETURN 
END 


t  303 
1  ~;3i 


SU0ROUTINC  sues 

BO  S 00 I  J*fSPPl #fSG«i 

f£U  JP  iJ } » f  Tgua  J  ( J  5  i  TABS  J  ?  J  }  «fXP  f  TAUOJ  C  J  i  /C1WSTAU2M  f  J!  I 

orrur^j 

0htS 


sot 

toe 

sas 
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SUmQi'TlUE  SU94 
1010  DC  1058 

I Qlt  ir(HFt!MD<  mtl0?2*l0!2»10!4 

105  2  SKFIINt  m*AKFI  < !  I  J^CTf***85CFS  (II  5  *£XP  <  ~Ci?  *  |  (  II 1/CT  pG&mFt  US  11 

1  Ol  3  <sO  TO  1018 

5014  TEMPI »0,0 

1013  00  1016  J*1*I5S*1 

1016  TEMPlsTEWPl  +BETA!JU1*J>  *  SIUSQTCJJ 

1017  SKFUNl  | !  } *EXP  C-TEMPi -CKFI  <  I  1  1/CTP**DKF!  1 1  1  » )*AKF I  1 1 1  l^CTP#* 
18KFH  1  !  }/<CONl*TEKPl»#lB£TAl <  f !  ) 

1018  CONTINUE 

1019  PETUftN 
END 


SI© 

-i 

ISO 

£ 

mb 

3 

Sl<§ 

m 

§ 

Si© 

& 

SI© 

f 

sid 

0 

sis 

9 

si© 

iJs 

sss 

n 

SIS 

S0 

SIO 

13 

037 


SUBROUTINE  3 US 5  $H 

1800  CALL  SSWTCHt  3  .KOOOrX)  Mi 

go  rgcioot 4soie>«KeoDFs  Sji 

1001  SFU^fW-fOi^nOOZ*  J002«!018  ill 

1002  white  isuoo^MSuyoore  jncepu?  jj  *i:psji 

i #XLAMJ< J?»CVUCJ1* J»? »13S*£ J  SI I 

1003  FORMAT  nP££l9«7ji  SIS 

1008  WHITE  1 2*  1009  HOF  IOC  m«CKP!  (II  )  k'KI  m  )*8KFI  fwu  1 1  *  SMfHiti 

HI  5t5St1un«CH2HUJ«!f"l*lSR«l  I  SH 

1009  P0HM<*T'tP7E15#7)  Si I 

1010  DO  1011  ll*ltISH«l  fit 

1011  WHITE  (2»1003HCQIJUI*J)»J*ISCPl»ISS«n  11! 

1017  WHITE  «2*  1003)  (BI  J*MU*  J«1*M«1)  Sj! 

1010  RETURN  SI! 

END  311 


4 


fit 

# 

7 

S 

# 

I© 

u 

If 

13 

14 

15 
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subroutine  sus* 

IF <  INDSUm  i  ' SO*  1 150*1000 
E 00©  00  1140  J«I«IS£*1 

ioio  ixxxeMSuswi  ji 
10 to  TT 1*0*0 
1030  TT£»0*0 
1040  ??3**0#0 

lO&w  00  1110  L*1*!XXX*1 

EC'SO  XX2*$G  JL. <  J  *L  ) *£XP  1  -CE  JL.P fJ«LJ/AITt211 

1 0T0  TT2*t V2+XX2 

1080  XXI*XX2*CEJLP( J*L> 

1090  TT1-TT1+XX1 

HOC  XX3*XX1  *CE  JLP<  J»t. ) 

1110  TT3**TT3+  XX3 

1120  EECHJC J>»TT!/TT2 

1130  £ECNUJ<JJsAL0GCT72/SGUL<J*  1  J  ) 

1140  E£CCPJ<J»a  (  TT2«TT3  -TTl *TT1 )  /((AIT (Z)  *  TT2 )  »*2  ) 
1150  RETURN 

END 


SUBROUTINE  TESTl 
ISO!  fFURKIND-l  >1060*1000*1060 
1000  tF(  I  VST  1 1 070* ' 010*1 070 
1010  00  1020  I»1 * ISR«! 

1 020  CHttTU  )«CHl 1  <  n 
IY3T*l 

1030  I F  (  YPRE V— A i  T  ( *  5 110*0*1060* 1 0*0 
1040  CALL  WRITE  k 

EXTEST* 1 .E-7*EXT£4J  <9  > 

NOSAYsO 
1060  RETURN 

*070  IF'AITUJ-YPREV)  1530.1060*1330 
1S30  DO  15*0  JiM«lSR*l 
15*0  CHI IT< It  *CHI 1(1) 

1550  DO  I 5"*Q  3*1*20*1 

1560  A*! ( J)=A*T(J} 

1570  A3! ( J )*A3T ( J 5 
1580  DO  1600  J-l *3* 1 
1590  A2 1 ( J 1 “A2T ( J ) 

1600  A1 1 < J 1 *A 1 T  C J 1 

1601  A2 HAS  *A2T (* 1 

ir(EXTRA3( 10)  .EQ.  0,0)  GO  TO  1620 

IF (  (EXTRAJt9)“-Al  1(1))  ,GT.  Ei’TEST)  GO  TO  1621 

EXTRA J  <  9 ) *EXTRA J 19) +EXTRA J ( 10) 

EXTEST* 1 .E-7*EXTRAJ  <9 ) 

1620  CALL  WRITE  2 

1621  iF(ABS(YSAV£-Al I f 1 » )  *LE«  EXTEST)  GO  TO  1622 
IFtNOSAY.NE.05  GO  TO  167* 

1622  NOSAY.O 

1630  1DELXC-10ELXG*  1 

I6AC  1FCIDELXC-IEXTU  >1670.1650*1650 

1650  IF (OELY— EXTRA J (15)1 1651 *1670*1670 

1651  DCLY*DELY*EXTRAJU2) 

1660  IDELXC*0 

1670  CONTINUE 

IF (EXTRA J( i 0 ) *EQ. 0. )G0  TO  1676 

IF((A1 1 (1 )+DELY-EXTRAJ(9) )  ,G£.  ( -EXTEST  ) )  GO  TO  1675 
1676  CONTINUE 
1673  RETURN 
167*  DELY-DYSAV 

nosay-o 
GO  TO  1670 
1675  NOSAYol 

YSAVE-AJ 1(1) 

DYSAVoDELY 

DELY«EXTRA J(9)-At  1(1) 

GO  TO  1673 
END 


SIS  14 
Si  3  12 
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